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

CoinOslFactorization2.cpp

/* $Id: CoinOslFactorization2.cpp 1207 2009-10-10 10:37:10Z forrest $ */
/* Copyright (C) 1987, 2009, International Business Machines
   Corporation and others.  All Rights Reserved. */
/*
  CLP_OSL - if defined use osl
  0 - don't unroll 2 and 3 - don't use in Gomory
  1 - don't unroll - do use in Gomory
  2 - unroll - don't use in Gomory
  3 - unroll and use in Gomory
*/
#include "CoinOslC.h"

#ifndef NDEBUG
extern int ets_count;
extern int ets_check;
#endif
#define COIN_REGISTER register
#define COIN_REGISTER2
#define COIN_REGISTER3 register
#ifdef COIN_USE_RESTRICT
# define COIN_RESTRICT2 __restrict
#else
# define COIN_RESTRICT2
#endif
static int c_ekkshfpo_scan2zero(COIN_REGISTER const EKKfactinfo * COIN_RESTRICT2 fact,const int * COIN_RESTRICT mpermu,
                   double *COIN_RESTRICT worki, double *COIN_RESTRICT worko, int * COIN_RESTRICT mptr)
{
  
  /* Local variables */
  int irow;
  double tolerance = fact->zeroTolerance;
  int nin=fact->nrow;
  int * COIN_RESTRICT mptrX=mptr;
  if ((nin&1)!=0) {
    irow=1;
    if (fact->packedMode) {
      int irow0= *mpermu;
      double dval;
      assert (irow0>=1&&irow0<=nin);
      mpermu++;
      dval=worki[irow0];
      if (NOT_ZERO(dval)) {
      worki[irow0]=0.0;
      if (fabs(dval) >= tolerance) {
        *(worko++)=dval;
        *(mptrX++) = 0;
      }
      }
    } else {
      int irow0= *mpermu;
      double dval;
      assert (irow0>=1&&irow0<=nin);
      mpermu++;
      dval=worki[irow0];
      if (NOT_ZERO(dval)) {
      worki[irow0]=0.0;
      if (fabs(dval) >= tolerance) {
        *worko=dval;
        *(mptrX++) = 0;
      }
      }
      worko++;
    }
  } else {
    irow=0;
  }
  if (fact->packedMode) {
    for (; irow < nin; irow+=2) {
      int irow0,irow1;
      double dval0,dval1;
      irow0=mpermu[0];
      irow1=mpermu[1];
      assert (irow0>=1&&irow0<=nin);
      assert (irow1>=1&&irow1<=nin);
      dval0=worki[irow0];
      dval1=worki[irow1];
      if (NOT_ZERO(dval0)) {
      worki[irow0]=0.0;
      if (fabs(dval0) >= tolerance) {
        *(worko++)=dval0;
        *(mptrX++) = irow+0;
      }
      }
      if (NOT_ZERO(dval1)) {
      worki[irow1]=0.0;
      if (fabs(dval1) >= tolerance) {
        *(worko++)=dval1;
        *(mptrX++) = irow+1;
      }
      }
      mpermu+=2;
    }
  } else {
    for (; irow < nin; irow+=2) {
      int irow0,irow1;
      double dval0,dval1;
      irow0=mpermu[0];
      irow1=mpermu[1];
      assert (irow0>=1&&irow0<=nin);
      assert (irow1>=1&&irow1<=nin);
      dval0=worki[irow0];
      dval1=worki[irow1];
      if (NOT_ZERO(dval0)) {
      worki[irow0]=0.0;
      if (fabs(dval0) >= tolerance) {
        worko[0]=dval0;
        *(mptrX++) = irow+0;
      }
      }
      if (NOT_ZERO(dval1)) {
      worki[irow1]=0.0;
      if (fabs(dval1) >= tolerance) {
        worko[1]=dval1;
        *(mptrX++) = irow+1;
      }
      }
      mpermu+=2;
      worko+=2;
    }
  }
  return static_cast<int>(mptrX-mptr);
}
/*
 * c_ekkshfpi_list executes the following loop:
 *
 * for (k=nincol, i=1; k; k--, i++) {
 *   int ipt = mptr[i];
 *   int irow = mpermu[ipt]; 
 *   worko[mpermu[irow]] = worki[i];
 *   worki[i] = 0.0;
 * }
 */
static int c_ekkshfpi_list(const int *COIN_RESTRICT mpermu, 
                     double *COIN_RESTRICT worki, 
                     double *COIN_RESTRICT worko,
                     const int * COIN_RESTRICT mptr, int nincol,
                     int * lastNonZero)
{
  int i,k,irow0,irow1;
  int first=COIN_INT_MAX;
  int last=0;
  /* worko was zeroed out outside */
  k = nincol;
  i = 0;
  if ((k&1)!=0) {
    int ipt=mptr[i];
    irow0=mpermu[ipt];
    first = CoinMin(irow0,first);
    last = CoinMax(irow0,last);
    i++;
    worko[irow0]=*worki;
    *worki++=0.0;
  }
  k=k>>1;
  for (; k; k--) {
    int ipt0 = mptr[i];
    int ipt1 = mptr[i+1];
    irow0 = mpermu[ipt0];
    irow1 = mpermu[ipt1];
    i+=2;
    first = CoinMin(irow0,first);
    last = CoinMax(irow0,last);
    first = CoinMin(irow1,first);
    last = CoinMax(irow1,last);
    worko[irow0] = worki[0];
    worko[irow1] = worki[1];
    worki[0]=0.0;
    worki[1]=0.0;
    worki+=2;
  }
  *lastNonZero=last;
  return first;
}
/*
 * c_ekkshfpi_list2 executes the following loop:
 *
 * for (k=nincol, i=1; k; k--, i++) {
 *   int ipt = mptr[i];
 *   int irow = mpermu[ipt]; 
 *   worko[mpermu[irow]] = worki[ipt];
 *   worki[ipt] = 0.0;
 * }
 */
static int c_ekkshfpi_list2(const int *COIN_RESTRICT mpermu, double *COIN_RESTRICT worki, double *COIN_RESTRICT worko,
                      const int * COIN_RESTRICT mptr, int nincol,
                     int * lastNonZero)
{
#if 1
  int i,k,irow0,irow1;
  int first=COIN_INT_MAX;
  int last=0;
  /* worko was zeroed out outside */
  k = nincol;
  i = 0;
  if ((k&1)!=0) {
    int ipt=mptr[i];
    irow0=mpermu[ipt];
    first = CoinMin(irow0,first);
    last = CoinMax(irow0,last);
    i++;
    worko[irow0]=worki[ipt];
    worki[ipt]=0.0;
  }
  k=k>>1;
  for (; k; k--) {
    int ipt0 = mptr[i];
    int ipt1 = mptr[i+1];
    irow0 = mpermu[ipt0];
    irow1 = mpermu[ipt1];
    i+=2;
    first = CoinMin(irow0,first);
    last = CoinMax(irow0,last);
    first = CoinMin(irow1,first);
    last = CoinMax(irow1,last);
    worko[irow0] = worki[ipt0];
    worko[irow1] = worki[ipt1];
    worki[ipt0]=0.0;
    worki[ipt1]=0.0;
  }
#else
  int first=COIN_INT_MAX;
  int last=0;
  /* worko was zeroed out outside */
  for (int i=0; i<nincol; i++) {
    int ipt = mptr[i];
    int irow = mpermu[ipt];
    first = CoinMin(irow,first);
    last = CoinMax(irow,last);
    worko[irow] = worki[ipt];
    worki[ipt]=0.0;
  }
#endif
  *lastNonZero=last;
  return first;
}
/*
 * c_ekkshfpi_list3 executes the following loop:
 *
 * for (k=nincol, i=1; k; k--, i++) {
 *   int ipt  = mptr[i];
 *   int irow = mpermu[ipt];
 *   worko[irow] = worki[i];
 *   worki[i] = 0.0;
 *   mptr[i] = mpermu[ipt];
 * }
 */
static void c_ekkshfpi_list3(const int *COIN_RESTRICT mpermu,
                double *COIN_RESTRICT worki, double *COIN_RESTRICT worko,
                int * COIN_RESTRICT mptr, int nincol)
{
  int i,k,irow0,irow1;
  /* worko was zeroed out outside */
  k = nincol;
  i = 0;
  if ((k&1)!=0) {
    int ipt=mptr[i];
    irow0=mpermu[ipt];
    mptr[i] = irow0;
    i++;
    worko[irow0]=*worki;
    *worki++=0.0;
  }
  k=k>>1;
  for (; k; k--) {
    int ipt0 = mptr[i];
    int ipt1 = mptr[i+1];
    irow0 = mpermu[ipt0];
    irow1 = mpermu[ipt1];
    mptr[i] = irow0;
    mptr[i+1] = irow1;
    i+=2;
    worko[irow0] = worki[0];
    worko[irow1] = worki[1];
    worki[0]=0.0;
    worki[1]=0.0;
    worki+=2;
  }
}
static int c_ekkscmv(COIN_REGISTER const EKKfactinfo * COIN_RESTRICT2 fact,int n, double *COIN_RESTRICT dwork, int *COIN_RESTRICT mptr, 
               double *COIN_RESTRICT dwork2)
{
  double tolerance = fact->zeroTolerance;
  int irow;
  const int * COIN_RESTRICT mptrsave = mptr;
  double * COIN_RESTRICT dwhere = dwork+1;
  if ((n&1)!=0) {
    if (NOT_ZERO(*dwhere)) {
      if (fabs(*dwhere) >= tolerance) {
      *++dwork2 = *dwhere;
      *++mptr = SHIFT_INDEX(1);
      } else {
      *dwhere = 0.0;
      }
    }
    dwhere++;
    irow=2;
  } else {
    irow=1;
  }
  for (n=n>>1;n;n--) {
    int second = NOT_ZERO(*(dwhere+1));
    if (NOT_ZERO(*dwhere)) {
      if (fabs(*dwhere) >= tolerance) {
      *++dwork2 = *dwhere;
      *++mptr = SHIFT_INDEX(irow);
      } else {
      *dwhere = 0.0;
      }
    }
    if (second) {
      if (fabs(*(dwhere+1)) >= tolerance) {
      *++dwork2 = *(dwhere+1);
      *++mptr = SHIFT_INDEX(irow+1);
      } else {
      *(dwhere+1) = 0.0;
      }
    }
    dwhere+=2;
    irow+=2;
  }
  
  return static_cast<int>(mptr-mptrsave);
} /* c_ekkscmv */
double c_ekkputl(const EKKfactinfo * COIN_RESTRICT2 fact,
           const int *COIN_RESTRICT mpt2,
           double *COIN_RESTRICT dwork1,
           double del3, 
           int nincol, int nuspik)
{
  double * COIN_RESTRICT dwork3     = fact->xeeadr+fact->nnentu;
  int * COIN_RESTRICT hrowi   = fact->xeradr+fact->nnentu;
  int offset = fact->R_etas_start[fact->nR_etas+1];
  int *COIN_RESTRICT hrowiR = fact->R_etas_index+offset;
  double *COIN_RESTRICT dluval = fact->R_etas_element+offset;
  int i, j;
  
  /* dwork1 is r', the new R transform
   * dwork3 is the updated incoming column, alpha_p
   * del3 apparently has the pivot of the incoming column (???).
   * Here, we compute the p'th element of R^-1 alpha_p
   * (as described on p. 273), which is just a dot product.
   * I don't know why we subtract.
   */
  for (i = 1; i <= nuspik; ++i) {
    j = UNSHIFT_INDEX(hrowi[ i]);
    del3 -= dwork3[i] * dwork1[j];
  }
  
  /* here we finally copy the r' to where we want it, the end */
  /* also take into account that the p'th row of R^-1 is -(p'th row of R). */
  /* also zero out dwork1 as we go */
  for (i = 0; i < nincol; ++i) {
    j = mpt2[i];
    hrowiR[ - i ] = SHIFT_INDEX(j);
    dluval[ - i ] = -dwork1[j];
    dwork1[j] = 0.;
  }
  
  return del3;
} /* c_ekkputl */
/* making this static seems to slow code down!
   may be being inlined
*/
int c_ekkputl2( const EKKfactinfo * COIN_RESTRICT2 fact,
           double *COIN_RESTRICT dwork1,
           double *del3p, 
           int nuspik)
{
  double * COIN_RESTRICT dwork3     = fact->xeeadr+fact->nnentu;
  int * COIN_RESTRICT hrowi   = fact->xeradr+fact->nnentu;
  int offset = fact->R_etas_start[fact->nR_etas+1];
  int *COIN_RESTRICT hrowiR = fact->R_etas_index+offset;
  double *COIN_RESTRICT dluval = fact->R_etas_element+offset;
  int i, j;
#if 0
  int nincol=c_ekksczr(fact,fact->nrow,dwork1,hrowiR);
#else
  int nrow=fact->nrow;
  const double tolerance = fact->zeroTolerance;
  int * COIN_RESTRICT mptrX=hrowiR;
  for (i = 1; i <= nrow; ++i) {
    if (dwork1[i] != 0.) {
      if (fabs(dwork1[i]) >= tolerance) {
      *(mptrX--) = SHIFT_INDEX(i);
      } else {
      dwork1[i] = 0.0;
      }
    }
  }
  int nincol=static_cast<int>(hrowiR-mptrX);
#endif
  double del3 = *del3p;
  /* dwork1 is r', the new R transform
   * dwork3 is the updated incoming column, alpha_p
   * del3 apparently has the pivot of the incoming column (???).
   * Here, we compute the p'th element of R^-1 alpha_p
   * (as described on p. 273), which is just a dot product.
   * I don't know why we subtract.
   */
  for (i = 1; i <= nuspik; ++i) {
    j = UNSHIFT_INDEX(hrowi[ i]);
    del3 -= dwork3[i] * dwork1[j];
  }
  
  /* here we finally copy the r' to where we want it, the end */
  /* also take into account that the p'th row of R^-1 is -(p'th row of R). */
  /* also zero out dwork1 as we go */
  for (i = 0; i < nincol; ++i) {
    j = UNSHIFT_INDEX(hrowiR[-i]);
    dluval[ - i ] = -dwork1[j];
    dwork1[j] = 0.;
  }
  
  *del3p = del3;
  return nincol;
} /* c_ekkputl */
static void c_ekkbtj4p_no_dense(const int nrow,const double * COIN_RESTRICT dluval, 
                        const int * COIN_RESTRICT hrowi,
                        const int * COIN_RESTRICT mcstrt,  
                        double * COIN_RESTRICT dwork1, int ndo,int jpiv)
{
  int i;
  double dv1;
  int iel;
  int irow;
  int i1,i2;
  
  /* count down to first nonzero */
  for (i=nrow;i >=1;i--) {
    if (dwork1[i]) {
      break;
    }
  }
  i--; /* as pivot is just 1.0 */
  if (i>ndo+jpiv) {
    i=ndo+jpiv;
  }
  mcstrt -= jpiv;
  i2=mcstrt[i+1];
  for (; i > jpiv; --i) {
    double dv1b=0.0;
    int nel;
    i1 = mcstrt[i];
    nel= i1-i2;
    dv1 = dwork1[i];
    iel=i2;
    if ((nel&1)!=0) {
      irow = hrowi[iel];
      dv1b = SHIFT_REF(dwork1, irow) * dluval[iel];
      iel++;
    }
    for ( ; iel < i1; iel+=2) {
      int irow = hrowi[iel];
      int irowb = hrowi[iel+1];
      dv1 += SHIFT_REF(dwork1, irow) * dluval[iel];
      dv1b += SHIFT_REF(dwork1, irowb) * dluval[iel+1];
    }
    i2=i1;
    dwork1[i] = dv1+dv1b;
  }
} /* c_ekkbtj4 */

static int c_ekkbtj4p_dense(const int nrow,const double * COIN_RESTRICT dluval,
                      const int * COIN_RESTRICT hrowi,
                      const int * COIN_RESTRICT mcstrt,  double * COIN_RESTRICT dwork1,
               int ndenuc,
               int ndo,int jpiv)
{
  int i;
  int i2;
  
  int last=ndo-ndenuc+1;
  double * COIN_RESTRICT densew = &dwork1[nrow-1];
  int nincol=0;
  const double * COIN_RESTRICT dlu1;
  dluval--;
  hrowi--;
  /* count down to first nonzero */
  for (i=nrow;i >=1;i--) {
    if (dwork1[i]) {
      break;
    }
  }
  if (i<ndo+jpiv) {
    int diff = ndo+jpiv-i;
    ndo -= diff;
    densew-=diff;
    nincol=diff;
  }
  i2=mcstrt[ndo+1];
  dlu1=&dluval[i2+1];
  for (i = ndo; i >last; i-=2) {
    int k;
    double dv1,dv2;
    const double * COIN_RESTRICT dlu2;
    dv1=densew[1];
    dlu2=dlu1+nincol;
    dv2=densew[0];
    for (k=0;k<nincol;k++) {
#ifdef DEBUG
      int kk=dlu1-dluval;
      int jj = (densew+(nincol-k+1))-dwork1;
      int ll=hrowi[k+kk];
      if (ll!=jj) abort();
#endif
      dv1 += densew[nincol-k+1]*dlu1[k];
      dv2 += densew[nincol-k+1]*dlu2[k];
    }
    densew[1]=dv1;
    dlu1=dlu2+nincol;
    dv2 += dv1*dlu1[0];
    dlu1++;
    nincol+=2;
    densew[0]=dv2;
    densew-=2;
  }
  return i;
} /* c_ekkbtj4 */

static void c_ekkbtj4p_after_dense(const double * COIN_RESTRICT dluval, 
                           const int * COIN_RESTRICT hrowi,
                           const int * COIN_RESTRICT mcstrt,  
                           double * COIN_RESTRICT dwork1, int i,int jpiv)
{
  int iel;
  mcstrt -= jpiv;
  i += jpiv;
  iel=mcstrt[i+1];
  for (; i > jpiv+1; i-=2) {
    int i1 = mcstrt[i];
    double dv1 = dwork1[i];
    double dv2;
    for (; iel < i1; iel++) {
      int irow = hrowi[iel];
      dv1 += SHIFT_REF(dwork1, irow) * dluval[iel];
    }
    i1 = mcstrt[i-1];
    dv2 = dwork1[i-1];
    dwork1[i] = dv1;
    for (; iel < i1; iel++) {
      int irow = hrowi[iel];
      dv2 += SHIFT_REF(dwork1, irow) * dluval[iel];
    }
    dwork1[i-1] = dv2;
  }
  if (i>jpiv) {
    int i1 = mcstrt[i];
    double dv1 = dwork1[i];
    for (; iel < i1; iel++) {
      int irow = hrowi[iel];
      dv1 += SHIFT_REF(dwork1, irow) * dluval[iel];
    }
    dwork1[i] = dv1;
  }
}

static void c_ekkbtj4p(COIN_REGISTER const EKKfactinfo * COIN_RESTRICT2 fact,
            double * COIN_RESTRICT dwork1)
{
  int lstart=fact->lstart;
  const int * COIN_RESTRICT hpivco  = fact->kcpadr;
  const double * COIN_RESTRICT dluval     = fact->xeeadr+1;
  const int * COIN_RESTRICT hrowi   = fact->xeradr+1;
  const int * COIN_RESTRICT mcstrt  = fact->xcsadr+lstart-1;
  int jpiv=hpivco[lstart]-1;
  int ndo=fact->xnetalval;
  /*     see if dense enough to unroll */
  if (fact->ndenuc<5) {
    c_ekkbtj4p_no_dense(fact->nrow,dluval,hrowi,mcstrt,dwork1,ndo,jpiv);
  } else {
    int i = c_ekkbtj4p_dense(fact->nrow,dluval,hrowi,mcstrt,dwork1,
                     fact->ndenuc, ndo,jpiv);
    c_ekkbtj4p_after_dense(dluval,hrowi,mcstrt,dwork1,i,jpiv);
  }
} /* c_ekkbtj4p */

static int c_ekkbtj4_sparse(COIN_REGISTER2 const EKKfactinfo * COIN_RESTRICT2 fact,
                    double * COIN_RESTRICT dwork1, 
                    int * COIN_RESTRICT mpt,    /* C style */
                    double * COIN_RESTRICT dworko,
                    int nincol, int * COIN_RESTRICT spare)
{
  const int nrow        = fact->nrow;
  const int * COIN_RESTRICT hcoli   = fact->xecadr;
  const int * COIN_RESTRICT mrstrt  = fact->xrsadr+nrow;
  char * COIN_RESTRICT nonzero = fact->nonzero;
  const int * COIN_RESTRICT hpivro  = fact->krpadr;
  const double * COIN_RESTRICT de2val     = fact->xe2adr-1;
  double tolerance = fact->zeroTolerance;
  double dv;
  int iel;
  
  int k,nStack,kx;
  int nList=0;
  int * COIN_RESTRICT list = spare;
  int * COIN_RESTRICT stack = spare+nrow;
  int * COIN_RESTRICT next = stack+nrow;
  int iPivot,kPivot;
  int iput,nput=0,kput=nrow;
  int j;
  int firstDoRow=fact->firstDoRow;
  
  for (k=0;k<nincol;k++) {
    nStack=1;
    iPivot=mpt[k];
    if (nonzero[iPivot]!=1&&iPivot>=firstDoRow) {
      stack[0]=iPivot;
      next[0]=mrstrt[iPivot];
      while (nStack) {
      /* take off stack */
      kPivot=stack[--nStack];
      if (nonzero[kPivot]!=1&&kPivot>=firstDoRow) {
        j=next[nStack];
        if (j==mrstrt[kPivot+1]) {
          /* finished so mark */
          list[nList++]=kPivot;
          nonzero[kPivot]=1;
        } else {
          kPivot=hcoli[j];
          /* put back on stack */
          next[nStack++] ++;
          if (!nonzero[kPivot]) {
            /* and new one */
            stack[nStack]=kPivot;
            nonzero[kPivot]=2;
            next[nStack++]=mrstrt[kPivot];
          }
        }
      } else if (kPivot<firstDoRow) {
        list[--kput]=kPivot;
        nonzero[kPivot]=1;
      }
      }
    } else if (nonzero[iPivot]!=1) {
      /* nothing to do (except check size at end) */
      list[--kput]=iPivot;
      nonzero[iPivot]=1;
    }
  }
  if (fact->packedMode) {
    dworko++;
    for (k=nList-1;k>=0;k--) {
      double dv;
      iPivot = list[k];
      dv = dwork1[iPivot];
      dwork1[iPivot]=0.0;
      nonzero[iPivot]=0;
      if (fabs(dv) > tolerance) {
      iput=hpivro[iPivot];
      kx=mrstrt[iPivot];
      dworko[nput]=dv;
      for (iel = kx; iel < mrstrt[iPivot+1]; iel++) {
        double dval;
        int irow = hcoli[iel];
        dval=de2val[iel];
        dwork1[irow] += dv*dval;
      }
      mpt[nput++]=iput-1;
      } else {
      dwork1[iPivot]=0.0;     /* force to zero, not just near zero */
      }
    }
    /* check remainder */
    for (k=kput;k<nrow;k++) {
      iPivot = list[k];
      nonzero[iPivot]=0;
      dv = dwork1[iPivot];
      dwork1[iPivot]=0.0;
      iput=hpivro[iPivot];
      if (fabs(dv) > tolerance) {
      dworko[nput]=dv;
      mpt[nput++]=iput-1;
      }
    }
  } else {
    /* not packed */
    for (k=nList-1;k>=0;k--) {
      double dv;
      iPivot = list[k];
      dv = dwork1[iPivot];
      dwork1[iPivot]=0.0;
      nonzero[iPivot]=0;
      if (fabs(dv) > tolerance) {
      iput=hpivro[iPivot];
      kx=mrstrt[iPivot];
      dworko[iput]=dv;
      for (iel = kx; iel < mrstrt[iPivot+1]; iel++) {
        double dval;
        int irow = hcoli[iel];
        dval=de2val[iel];
        dwork1[irow] += dv*dval;
      }
      mpt[nput++]=iput-1;
      } else {
      dwork1[iPivot]=0.0;     /* force to zero, not just near zero */
      }
    }
    /* check remainder */
    for (k=kput;k<nrow;k++) {
      iPivot = list[k];
      nonzero[iPivot]=0;
      dv = dwork1[iPivot];
      dwork1[iPivot]=0.0;
      iput=hpivro[iPivot];
      if (fabs(dv) > tolerance) {
      dworko[iput]=dv;
      mpt[nput++]=iput-1;
      }
    }
  }
  
  return (nput);
} /* c_ekkbtj4 */

static void c_ekkbtjl(COIN_REGISTER2 const EKKfactinfo * COIN_RESTRICT2 fact,
                double * COIN_RESTRICT dwork1)
{
  int i, j, k, k1;
  int l1;
  const double * COIN_RESTRICT dluval = fact->R_etas_element;
  const int * COIN_RESTRICT hrowi = fact->R_etas_index;
  const int * COIN_RESTRICT mcstrt = fact->R_etas_start;
  const int * COIN_RESTRICT hpivco = fact->hpivcoR;
  int ndo=fact->nR_etas;
#ifndef UNROLL1
#define UNROLL1 4
#endif
#if UNROLL1>2
  int l2;
#endif
  int kn;
  double dv;
  int iel;
  int ipiv;
  int knext;
  
  knext = mcstrt[ndo + 1];
#if UNROLL1>2
  for (i = ndo; i > 0; --i) {
    k1 = knext;
    knext = mcstrt[i];
    ipiv = hpivco[i];
    dv = dwork1[ipiv];
    /*       fast floating */
    k = knext - k1;
    kn = k >> 2;
    iel = k1 + 1;
    if (dv != 0.) {
      l1 = (k & 1) != 0;
      l2 = (k & 2) != 0;
      for (j = 1; j <= kn; j++) {
      int irow0 = hrowi[iel + 0];
      int irow1 = hrowi[iel + 1];
      int irow2 = hrowi[iel + 2];
      int irow3 = hrowi[iel + 3];
      double dval0 = dv * dluval[iel + 0] + SHIFT_REF(dwork1, irow0);
      double dval1 = dv * dluval[iel + 1] + SHIFT_REF(dwork1, irow1);
      double dval2 = dv * dluval[iel + 2] + SHIFT_REF(dwork1, irow2);
      double dval3 = dv * dluval[iel + 3] + SHIFT_REF(dwork1, irow3);
      SHIFT_REF(dwork1, irow0) = dval0;
      SHIFT_REF(dwork1, irow1) = dval1;
      SHIFT_REF(dwork1, irow2) = dval2;
      SHIFT_REF(dwork1, irow3) = dval3;
      iel+=4;
      }
      if (l1) {
      int irow0 = hrowi[iel];
      SHIFT_REF(dwork1, irow0) += dv* dluval[iel];
      ++iel;
      }
      if (l2) {
      int irow0 = hrowi[iel + 0];
      int irow1 = hrowi[iel + 1];
      SHIFT_REF(dwork1, irow0) += dv* dluval[iel];
      SHIFT_REF(dwork1, irow1) += dv* dluval[iel+1];
      }
    }
  }
#else
  for (i = ndo; i > 0; --i) {
    k1 = knext;
    knext = mcstrt[i];
    ipiv = hpivco[i];
    dv = dwork1[ipiv];
    k = knext - k1;
    kn = k >> 1;
    iel = k1 + 1;
    if (dv != 0.) {
      l1 = (k & 1) != 0;
      for (j = 1; j <= kn; j++) {
      int irow0 = hrowi[iel + 0];
      int irow1 = hrowi[iel + 1];
      double dval0 = dv * dluval[iel + 0] + SHIFT_REF(dwork1, irow0);
      double dval1 = dv * dluval[iel + 1] + SHIFT_REF(dwork1, irow1);
      SHIFT_REF(dwork1, irow0) = dval0;
      SHIFT_REF(dwork1, irow1) = dval1;
      iel+=2;
      }
      if (l1) {
      int irow0 = hrowi[iel];
      SHIFT_REF(dwork1, irow0) += dv* dluval[iel];
      ++iel;
      }
    }
  }
#endif
} /* c_ekkbtjl */

static int c_ekkbtjl_sparse(COIN_REGISTER2 const EKKfactinfo * COIN_RESTRICT2 fact,
               double * COIN_RESTRICT dwork1, 
               int * COIN_RESTRICT mpt , int nincol)
{
  const double * COIN_RESTRICT dluval = fact->R_etas_element;
  const int * COIN_RESTRICT hrowi = fact->R_etas_index;
  const int * COIN_RESTRICT mcstrt = fact->R_etas_start;
  const int * COIN_RESTRICT hpivco = fact->hpivcoR;
  char * COIN_RESTRICT nonzero = fact->nonzero;
  int ndo=fact->nR_etas;
  int i, j, k1;
  double dv;
  int ipiv;
  int irow0, irow1;
  int knext;
  int number=nincol;
  
  /*     ------------------------------------------- */
  /* adjust back */
  hrowi++;
  dluval++;
  
  /*         DO ANY ROW TRANSFORMATIONS */
  
  /* Function Body */
  knext = mcstrt[ndo + 1];
  for (i = ndo; i > 0; --i) {
    k1 = knext;
    knext = mcstrt[i];
    ipiv = hpivco[i];
    dv = dwork1[ipiv];
    if (dv) {
      for (j = k1; j <knext-1; j+=2) {
      irow0 = hrowi[j];
      irow1 = hrowi[j+1];
      SHIFT_REF(dwork1, irow0) += dv * dluval[j];
      SHIFT_REF(dwork1, irow1) += dv * dluval[j+1];
      if (!nonzero[irow0]) {
        nonzero[irow0]=1;
        mpt[++number]=UNSHIFT_INDEX(irow0);
      }
      if (!nonzero[irow1]) {
        nonzero[irow1]=1;
        mpt[++number]=UNSHIFT_INDEX(irow1);
      }
      }
      if (j<knext) {
      irow0 = hrowi[j];
      SHIFT_REF(dwork1, irow0) += dv * dluval[j];
      if (!nonzero[irow0]) {
        nonzero[irow0]=1;
        mpt[++number]=UNSHIFT_INDEX(irow0);
      }
      }
    }
  }
  return (number);
} /* c_ekkbtjl */



static void c_ekkbtju_dense(const int nrow,
                      const double * COIN_RESTRICT dluval,
                      const int * COIN_RESTRICT hrowi,
                      const int * COIN_RESTRICT mcstrt, 
                      int * COIN_RESTRICT hpivco, 
                      double * COIN_RESTRICT dwork1,
                      int * COIN_RESTRICT start,int last,int offset,
                      double * COIN_RESTRICT densew)
{
  /* Local variables */
  int ipiv1,ipiv2;
  int save=hpivco[last];
  
  hpivco[last]=nrow+1;
  
  ipiv1=*start;
  ipiv2=hpivco[ipiv1];
  while(ipiv2<last) {
    int iel,k;
    const int   kx1     = mcstrt[ipiv1];
    const int   kx2  = mcstrt[ipiv2];
    const int   nel1 = hrowi[kx1-1];
    const int   nel2 = hrowi[kx2-1];
    const double dpiv1 = dluval[kx1-1];
    const double dpiv2 = dluval[kx2-1];
    const int   n1      = offset+ipiv1; /* number in dense part */
    const int   nsparse1=nel1-n1;
    const int   nsparse2=nel2-n1-(ipiv2-ipiv1);
    const int   k1 = kx1+nsparse1;
    const int   k2 = kx2+nsparse2;
    const double *dlu1 = &dluval[k1];
    const double *dlu2 = &dluval[k2];
    
    double dv1 = dwork1[ipiv1];
    double dv2 = dwork1[ipiv2];
    
    for (iel = kx1; iel < k1; ++iel) {
      dv1 -= SHIFT_REF(dwork1, hrowi[iel]) * dluval[iel];
    }
    for (iel = kx2; iel < k2; ++iel) {
      dv2 -= SHIFT_REF(dwork1, hrowi[iel]) * dluval[iel];
    }
    for (k=0;k<n1;k++) {
      dv1 -= dlu1[k] * densew[k];
      dv2 -= dlu2[k] * densew[k];
    }
    dv1 *= dpiv1;
    dv2 -= dlu2[n1] * dv1;
    dwork1[ipiv1] = dv1;
    dwork1[ipiv2] = dv2*dpiv2;
    ipiv1 = hpivco[ipiv2];
    ipiv2 = hpivco[ipiv1];
  }
  hpivco[last]=save;
  
  *start=ipiv1;
  return;
}
/* about 8-10% of execution time is spent in this routine */
static int c_ekkbtju_aux(const double * COIN_RESTRICT dluval, 
                   const int * COIN_RESTRICT hrowi,
                   const int * COIN_RESTRICT mcstrt, 
                   const int * COIN_RESTRICT hpivco,
                   double * COIN_RESTRICT dwork1,
                   int ipiv, int loop_end)
{
#define UNROLL2 2
#ifndef UNROLL2
#if CLP_OSL==2||CLP_OSL==3
#define UNROLL2 2
#else
#define UNROLL2 1
#endif
#endif
  while (ipiv<=loop_end) {
    int kx = mcstrt[ipiv];
    const int nel = hrowi[kx-1];
#if UNROLL2<2
    const int kxe = kx + nel;
#endif
    
    double dv = dwork1[ipiv]; /* rhs */
#if UNROLL2>1
    const int * hrowi2=hrowi+kx;
    const int * hrowi2end=hrowi2+nel;
    const double * dluval2=dluval+kx;
#else
    int iel;
#endif
    const double dpiv = dluval[kx-1];     /* inverse of pivot */
    
    
    /* subtract terms whose unknowns have been solved for */
    
    /* a significant proportion of these loops may not modify dv at all.
     * However, it seems to be just as expensive to check if the loop
     * would modify dv as it is to just do it.
     * The only difference would be that dluval wouldn't be referenced
     * for those loops, would might save some cache paging,
     * but unfortunately the code generated to search for zeros (on AIX)
     * is *worse* than code that just multiplies by dval.
     */
#if UNROLL2<2
    for (iel = kx; iel < kxe; ++iel) {
      const int irow = hrowi[iel];
      const double dval=dluval[iel];
      dv -= SHIFT_REF(dwork1, irow) * dval;
    }
    
    dwork1[ipiv] = dv * dpiv; /* divide by the pivot */
#else
    if ((nel&1)!=0) {
      int irow = *hrowi2;
      double dval=*dluval2;
      dv -= SHIFT_REF(dwork1, irow) * dval;
      hrowi2++;
      dluval2++;
    }
    for (; hrowi2 < hrowi2end; hrowi2 +=2,dluval2 +=2) {
      int irow0 = hrowi2[0];
      int irow1 = hrowi2[1];
      double dval0=dluval2[0];
      double dval1=dluval2[1];
      double d0=SHIFT_REF(dwork1, irow0);
      double d1=SHIFT_REF(dwork1, irow1);
      dv -= d0 * dval0;
      dv -= d1 * dval1;
    }
    dwork1[ipiv] = dv * dpiv; /* divide by the pivot */
#endif
    
    ipiv=hpivco[ipiv];
  }
  
  return (ipiv);
}

/*
 * We are given the upper diagonal matrix U from the LU factorization
 * and a rhs dwork1.
 * This solves the system U x = dwork1
 * by back substitution, overwriting dwork1 with the solution x.
 *
 * It does this in textbook style by solving the equations "bottom" up,
 * so for each equation one new unknown is solved for by subtracting
 * from the rhs the sum of the terms whose unknowns have been solved for,
 * then dividing by the coefficient of the new unknown.
 *
 * Since we update the U matrix using F-T, the order of the columns
 * changes slightly each iteration.  Initially, hpivco[i] == i+1,
 * and each iteration (generally) introduces one element where this
 * is no longer true.  However, because we periodically refactorize,
 * it is much more common for hpivco[i] == i+1 than not.
 *
 * The one quirk is that value referred to as the pivot is actually
 * the reciprocal of the pivot, to avoid a division.
 *
 * Solving in this fashion is inappropriate if there are frequently
 * cases where all unknowns in an equation have value zero.
 * This seems to happen frequently if the sparsity of the rhs is, say, 10%.
 */
static void c_ekkbtju(COIN_REGISTER2 const EKKfactinfo * COIN_RESTRICT2 fact,      
           double * COIN_RESTRICT dwork1,
           int ipiv)
{
  const int nrow  = fact->nrow;
  double * COIN_RESTRICT dluval     = fact->xeeadr;
  int * COIN_RESTRICT hrowi   = fact->xeradr;
  int * COIN_RESTRICT mcstrt  = fact->xcsadr;
  int * COIN_RESTRICT hpivco_new = fact->kcpadr+1;
  int ndenuc=fact->ndenuc;
  int first_dense = fact->first_dense;
  int last_dense = fact->last_dense;
  
  const int has_dense = (first_dense<last_dense &&
                   mcstrt[ipiv]<=mcstrt[last_dense]);
  
  /* Parameter adjustments */
  /* dluval and hrowi were NOT decremented here.
     I believe that they are used as C-style arrays below.
     At this point, I am going to convert them from Fortran- to C-style
     here by incrementing them; at some later time, I will convert their
     uses in this file to Fortran-style.
  */
  dluval++;
  hrowi++;
  
  if (has_dense) 
    ipiv = c_ekkbtju_aux(dluval, hrowi, mcstrt, hpivco_new, dwork1, ipiv,
                   first_dense - 1);
  
  if (has_dense) {
    int n=0;
    int firstDense = nrow-ndenuc+1;
    double *densew = &dwork1[firstDense];
    
    /* check first dense to see where in triangle it is */
    int last=first_dense;
    int j=mcstrt[last]-1;
    int k1=j;
    int k2=j+hrowi[j];
    
    for (j=k2;j>k1;j--) {
      int irow=UNSHIFT_INDEX(hrowi[j]);
      if (irow<firstDense) {
      break;
      } else {
#ifdef DEBUG
      if (irow!=last-1) {
        abort();
      }
#endif
      last=irow;
      n++;
      }
    }
    c_ekkbtju_dense(nrow,dluval,hrowi,mcstrt,const_cast<int *> (hpivco_new),
              dwork1,&ipiv,last_dense, n - first_dense, densew);
  }
  
  (void) c_ekkbtju_aux(dluval, hrowi, mcstrt, hpivco_new, dwork1, ipiv, nrow);
} /* c_ekkbtju */


/*
 * mpt / *nincolp contain the indices of nonzeros in dwork1.
 * nonzero contains the same information as a byte-mask.
 *
 * currently, erase_nonzero is true iff this is called from c_ekketsj.
 */
static int c_ekkbtju_sparse(COIN_REGISTER2 const EKKfactinfo * COIN_RESTRICT2 fact,
               double * COIN_RESTRICT dwork1,
               int * COIN_RESTRICT mpt, int nincol,
               int * COIN_RESTRICT spare)
{
  const double * COIN_RESTRICT dluval     = fact->xeeadr+1;
  const int * COIN_RESTRICT mcstrt  = fact->xcsadr;
  char * COIN_RESTRICT nonzero = fact->nonzero;
  const int * COIN_RESTRICT hcoli   = fact->xecadr;
  const int * COIN_RESTRICT mrstrt  = fact->xrsadr;
  const int * COIN_RESTRICT hinrow  = fact->xrnadr;
  const double * COIN_RESTRICT de2val     = fact->xe2adr-1;
  int i;
  int iPivot;
  int nList=0;
  int nStack,k,kx;
  const int nrow=fact->nrow;
  const double tolerance = fact->zeroTolerance;
  int * COIN_RESTRICT list = spare;
  int * COIN_RESTRICT stack = spare+nrow;
  int * COIN_RESTRICT next = stack+nrow;
  /*
   * Examine all nonzero elements and determine which elements may be
   * nonzero in the result.
   * Any row in U that contains terms that may have nonzero variable values
   * may produce a nonzero value.
   */
  for (k=0;k<nincol;k++) {
    nStack=1;
    iPivot=mpt[k];
    stack[0]=iPivot;
    next[0]=0;
    while (nStack) {
      int kPivot,ninrow,j;
      /* take off stack */
      kPivot=stack[--nStack];
      /*printf("nStack %d kPivot %d, ninrow %d, j %d, nList %d\n",
      nStack,kPivot,hinrow[kPivot],
      next[nStack],nList);*/
      if (nonzero[kPivot]!=1) {
      ninrow = hinrow[kPivot];
      j=next[nStack];
      if (j!=ninrow) {
        kx = mrstrt[kPivot];
        kPivot=hcoli[kx+j];
        /* put back on stack */
        next[nStack++] ++;
        if (!nonzero[kPivot]) {
          /* and new one */
          stack[nStack]=kPivot;
          nonzero[kPivot]=2;
          next[nStack++]=0;
        }
      } else {
        /* finished so mark */
        list[nList++]=kPivot;
        nonzero[kPivot]=1;
      }
      }
    }
  }
  
  i=nList-1;
  nList=0;
  for (;i>=0;i--) {
    double dpiv;
    double dv;
    iPivot = list[i];
    kx      = mcstrt[iPivot];
    dpiv = dluval[kx-1];
    dv      = dpiv * dwork1[iPivot];
    nonzero[iPivot] = 0;
    if (fabs(dv)>=tolerance) {
      int iel;
      int krx     = mrstrt[iPivot];
      int krxe    = krx+hinrow[iPivot];
      dwork1[iPivot]=dv;
      mpt[nList++]=iPivot;
      for (iel = krx; iel < krxe; iel++) {
      int irow0 = hcoli[iel];
      double dval=de2val[iel];
      dwork1[irow0] -= dv*dval;
      }
    } else {
      dwork1[iPivot]=0.0;
    }
  }
  
  return (nList);
} /* c_ekkbtjuRow */

/*
 * dpermu is supposed to be zeroed on entry to this routine.
 * It is used as a working buffer.
 * The input vector dwork1 is permuted into dpermu, operated on,
 * and the answer is permuted back into dwork1, zeroing dpermu in
 * the process.
 */
/*
 * nincol > 0 ==> mpt contains indices of non-zeros in dpermu
 *
 * first_nonzero contains index of first (last??)nonzero;
 * only used if nincol==0.
 *
 * dpermu contains permuted input; dwork1 is now zero
 */
int c_ekkbtrn(COIN_REGISTER3 const EKKfactinfo * COIN_RESTRICT2 fact,
          double * COIN_RESTRICT dwork1,
          int * COIN_RESTRICT mpt, int first_nonzero)
{
  double * COIN_RESTRICT dpermu = fact->kadrpm;
  const int * COIN_RESTRICT mpermu=fact->mpermu;
  const int * COIN_RESTRICT hpivco_new= fact->kcpadr+1;
  
  const int nrow  = fact->nrow;
  int i;
  int nincol;
  /* find the first non-zero input */
  int ipiv;
  if (first_nonzero) {
    ipiv = first_nonzero;
#if 1
    if (c_ekk_IsSet(fact->bitArray,ipiv)) {
      /* slack */
      int lastSlack = fact->lastSlack;
      int firstDo=hpivco_new[lastSlack];
      assert (dpermu[ipiv]);
      while (ipiv!=firstDo) {
      assert (c_ekk_IsSet(fact->bitArray,ipiv));
      if (dpermu[ipiv])
        dpermu[ipiv]=-dpermu[ipiv];
      ipiv=hpivco_new[ipiv];
      }
    }
#endif
  } else {
    int lastSlack = fact->numberSlacks;
    ipiv=hpivco_new[0];
    for (i=0;i<lastSlack;i++) {
      int next_piv = hpivco_new[ipiv];
      assert (c_ekk_IsSet(fact->bitArray,ipiv));
      if (dpermu[ipiv]) {
      break;
      } else {
      ipiv=next_piv;
      }
    }
    
    /* usually, there is a non-zero slack entry... */
    if (i==lastSlack) {
      /* but if there isn't... */
      for (;i<nrow;i++) {
      if (!dpermu[ipiv]) {
        ipiv=hpivco_new[ipiv];
      } else {
        break;
      }
      }
    } else {
      /* reverse signs for slacks */
      for (;i<lastSlack;i++) {
      assert (c_ekk_IsSet(fact->bitArray,ipiv));
      if (dpermu[ipiv])
        dpermu[ipiv]=-dpermu[ipiv];
      ipiv=hpivco_new[ipiv];
      }
      assert (!c_ekk_IsSet(fact->bitArray,ipiv)||ipiv>fact->nrow);
      
      /* this is presumably the first non-zero non slack */
      /*ipiv=firstDo;*/
    }
  }
  if (ipiv<=fact->nrow) {
    /* skipBtju is always (?) 0 first the first call,
     * ipiv tends to be >nrow for the second */
    
    /*       DO U */
    c_ekkbtju(fact,dpermu,
          ipiv); 
  }
  
  
  /*       DO ROW ETAS IN L */
  c_ekkbtjl(fact, dpermu); 
  c_ekkbtj4p(fact,dpermu);
  
  /* dwork1[mpermu] = dpermu; dpermu = 0; mpt = indices of non-zeros */
  nincol = 
    c_ekkshfpo_scan2zero(fact,&mpermu[1],dpermu,&dwork1[1],&mpt[1]);
  
  /* dpermu should be zero now */
#ifdef DEBUG
  for (i=1;i<=nrow ;i++ ) {
    if (dpermu[i]) {
      abort();
    } /* endif */
  } /* endfor */
#endif
  return (nincol);
} /* c_ekkbtrn */

static int c_ekkbtrn0_new(COIN_REGISTER3 const EKKfactinfo * COIN_RESTRICT2 fact,
                  double * COIN_RESTRICT dwork1,
                  int * COIN_RESTRICT mpt, int nincol, 
                     int * COIN_RESTRICT spare)
{
  double * COIN_RESTRICT dpermu = fact->kadrpm;
  const int * COIN_RESTRICT mpermu=fact->mpermu;
  const int * COIN_RESTRICT hpivro  = fact->krpadr;
  
  const int nrow  = fact->nrow;
  
  int i;
  char * nonzero=fact->nonzero;
  int doSparse=1;
  
  /* so:  dpermu must contain room for:
   * nrow doubles, followed by
   * nrow ints (mpermu), followed by
   * nrow ints (the inverse permutation), followed by
   * an unused area (?) of nrow ints, followed by
   * nrow chars (this non-zero array).
   *
   * and apparently the first nrow elements of nonzero are expected
   * to already be zero.
   */    
#ifdef DEBUG
  for (i=1;i<=nrow ;i++ ) {
    if (nonzero[i]) {
      abort();
    } /* endif */
  } /* endfor */
#endif
  /* now nonzero[i]==1 iff there is an entry for i in mpt */
  
  nincol=c_ekkbtju_sparse(fact, dpermu,
                  &mpt[1], nincol,
                  spare);
  
  /* the vector may have more nonzero elements now */
  /*       DO ROW ETAS IN L */
#define DENSE_THRESHOLD (nincol*10+100)
  if (DENSE_THRESHOLD>nrow) {
    doSparse=0;
    c_ekkbtjl(fact, dpermu); 
  } else {
    /* set nonzero */
    for(i=0;i<nincol;i++) {
      int j=mpt[i+1];
      nonzero[j]=1;
    }
    nincol =
      c_ekkbtjl_sparse(fact,
                 dpermu,
                 mpt, 
                 nincol);
    for(i=0;i<nincol;i++) {
      int j=mpt[i+1];
      nonzero[j]=0;
    }
    if (DENSE_THRESHOLD>nrow) {
      doSparse=0;
#ifdef DEBUG
      for (i=1;i<=nrow;i++) {
      if (nonzero[i]) {
        abort();
      }
      }
#endif
    }
  }
  if (!doSparse) {
    c_ekkbtj4p(fact,dpermu);
    /* dwork1[mpermu] = dpermu; dpermu = 0; mpt = indices of non-zeros */
    nincol = 
      c_ekkshfpo_scan2zero(fact,&mpermu[1],dpermu,&dwork1[1],&mpt[1]);
    
    /* dpermu should be zero now */
#ifdef DEBUG
    for (i=1;i<=nrow ;i++ ) {
      if (dpermu[i]) {
      abort();
      } /* endif */
    } /* endfor */
#endif
  } else {
    /* still sparse */
    if (fact->nnentl) {
      nincol =
      c_ekkbtj4_sparse(fact,
                   dpermu,
                   &mpt[1],
                   dwork1,
                   nincol,spare);
    } else {
      double tolerance=fact->zeroTolerance;
      int irow;
      int nput=0;
      if (fact->packedMode) {
      for (i = 0; i <nincol; i++) {
        int irow0;
        double dval;
        irow=mpt[i+1];
        dval=dpermu[irow];
        if (NOT_ZERO(dval)) {
          if (fabs(dval) >= tolerance) {
            irow0= hpivro[irow];
            dwork1[1+nput]=dval;
            mpt[1 + nput++]=irow0-1;
          }
          dpermu[irow]=0.0;
        }
      }
      } else {
      for (i = 0; i <nincol; i++) {
        int irow0;
        double dval;
        irow=mpt[i+1];
        dval=dpermu[irow];
        if (NOT_ZERO(dval)) {
          if (fabs(dval) >= tolerance) {
            irow0= hpivro[irow];
            dwork1[irow0]=dval;
            mpt[1 + nput++]=irow0-1;
          }
          dpermu[irow]=0.0;
        }
      }
      }
      nincol=nput;
    }
  }
  
  
  return (nincol);
} /* c_ekkbtrn */


/* returns c_ekkbtrn(fact, dwork1, mpt)
 *
 * but since mpt[1..nincol] contains the indices of non-zeros in dwork1,
 * we can do faster.
 */
static int c_ekkbtrn_mpt(COIN_REGISTER3 const EKKfactinfo * COIN_RESTRICT2 fact,
            double * COIN_RESTRICT dwork1,
            int * COIN_RESTRICT mpt, int nincol,int * COIN_RESTRICT spare)
{
  double * COIN_RESTRICT dpermu = fact->kadrpm;
  const int nrow  = fact->nrow;
  
  const int * COIN_RESTRICT mpermu=fact->mpermu;
  /*const int *mrstrt   = fact->xrsadr;*/
  
#ifdef DEBUG
  int i;
  memset(spare,'A',3*nrow*sizeof(int));
  {
    
    for (i=1;i<=nrow;i++) {
      if (dpermu[i]) {
      abort();
      }
    }
  } 
#endif
  
  
  int i;
#ifdef DEBUG
  for (i=1;i<=nrow;i++) {
    if (fact->nonzero[i]||dpermu[i]) {
      abort();
    }
  }
#endif
  assert (fact->if_sparse_update>0&&mpt&&fact->rows_ok) ;
  
  /* read the input vector from mpt/dwork1;
   * permute it into dpermu;
   * construct a nonzero mask in nonzero;
   * overwrite mpt with the permuted indices;
   * clear the dwork1 vector.
   */
  for (i=0;i<nincol;i++) {
    int irow=mpt[i+1];
    int jrow=mpermu[irow];
    dpermu[jrow]=dwork1[irow];
    /*nonzero[jrow-1]=1; this is done in btrn0 */
    mpt[i+1]=jrow;
    dwork1[irow]=0.0;
  }
  
  if (DENSE_THRESHOLD<nrow) {
    nincol = c_ekkbtrn0_new(fact, dwork1, mpt, nincol,spare);
  } else {
    nincol = c_ekkbtrn(fact, dwork1, mpt, 0);
  }
#ifdef DEBUG
  {
    
    for (i=1;i<=nrow;i++) {
      if (dpermu[i]) {
      abort();
      }
    }
    if (fact->if_sparse_update>0) {
      for (i=1;i<=nrow;i++) {
      if (fact->nonzero[i]) {
        abort();
      }
      }
    }
  } 
#endif
  return nincol;
}

/* returns c_ekkbtrn(fact, dwork1, mpt)
 *
 * but since (dwork1[i]!=0) == (i==ipivrw),
 * we can do faster.
 */
int c_ekkbtrn_ipivrw(COIN_REGISTER3 const EKKfactinfo * COIN_RESTRICT2 fact,
               double * COIN_RESTRICT dwork1,
               int * COIN_RESTRICT mpt, int ipivrw,int * COIN_RESTRICT spare)
{
  double * COIN_RESTRICT dpermu = fact->kadrpm;
  const int nrow  = fact->nrow;
  
  const int * COIN_RESTRICT mpermu=fact->mpermu;
  const double * COIN_RESTRICT dluval     = fact->xeeadr;
  const int * COIN_RESTRICT mrstrt  = fact->xrsadr;
  const int * COIN_RESTRICT hinrow  = fact->xrnadr;
  const int * COIN_RESTRICT hcoli   = fact->xecadr;
  const int * COIN_RESTRICT mcstrt  = fact->xcsadr;
  
  int nincol;
  
#ifdef DEBUG
  int i;
  for (i=1;i<=nrow ;i++ ) {
    if (dpermu[i]) {
      abort();
    } /* endif */
  } /* endfor */
#endif
  
  if (fact->if_sparse_update>0&&mpt&& fact->rows_ok) {
    mpt[1] = ipivrw;
    nincol = c_ekkbtrn_mpt(fact, dwork1, mpt, 1,spare);
  } else {
    int ipiv;
    int kpivrw = mpermu[ipivrw];
    dpermu[kpivrw]=dwork1[ipivrw];
    dwork1[ipivrw]=0.0;
    
    if (fact->rows_ok) {
      /* !fact->if_sparse_update
       * but we still have rowwise info,
       * so we may as well use it to do the slack row
       */
      int iipivrw=nrow+1;
      int itest = fact->nnentu+1;
      int k=mrstrt[kpivrw];
      int lastInRow= k+hinrow[kpivrw];
      double dpiv,dv;
      for (;k<lastInRow;k++) {
      int icol=hcoli[k];
      int start=mcstrt[icol];
      if (start<itest) {
        iipivrw=icol;
        itest=start;
      }
      }
      /* do missed pivot */
      itest=mcstrt[kpivrw];
      dpiv=dluval[itest];
      dv=dpermu[kpivrw];
      dv*=dpiv;
      dpermu[kpivrw]=dv;
      ipiv=iipivrw;
    } else {
      /* no luck - c_ekkbtju will slog through slacks (?) */
      ipiv=kpivrw;
    }
    /* nincol not read */
    /* not sparse */
    /* do slacks */
    if (ipiv<=fact->nrow) {
      if (c_ekk_IsSet(fact->bitArray,ipiv)) {
      const int * hpivco_new= fact->kcpadr+1;
      int lastSlack = fact->lastSlack;
      int firstDo=hpivco_new[lastSlack];
      /* slack */
      /* need pivot row of first nonslack */
      dpermu[ipiv]=-dpermu[ipiv];
#ifndef NDEBUG
      while (1) {
        assert (c_ekk_IsSet(fact->bitArray,ipiv));
        ipiv=hpivco_new[ipiv];
        if (ipiv>fact->nrow||ipiv==firstDo)
          break;
      }
      assert (!c_ekk_IsSet(fact->bitArray,ipiv)||ipiv>fact->nrow);
      assert (ipiv==firstDo);
#endif
      ipiv=firstDo;
      }
    }
    nincol = c_ekkbtrn(fact, dwork1, mpt, ipiv);
  }
  
  return nincol;
}
/*
 * Does work associated with eq. 3.7:
 *    r' = u' U^-1
 *
 * where u' (can't write the overbar) is the p'th row of U, without
 * the entry for column p.  (here, jpivrw is p).
 * We solve this as for btju.  We know
 *    r' U = u'
 *
 * so we solve from low index to hi, determining the next value u_i'
 * by doing the dot-product of r' and the i'th column of U (excluding
 * element i itself), subtracting that from u'_i, and dividing by
 * U_ii (we store the reciprocal, so here we multiply).
 *
 * Now, in principle dwork1 should be initialized to the p'th row of U.
 * Instead, it is initially zeroed and filled in as we go along.
 * Of the entries in u' that we reference during a dot product with
 * a column of U, either
 *    the entry is 0 by definition, since it is < p, or
 *    it has already been set by a previous iteration, or
 *    it is p.
 *
 * Because of this, we know that all elements < p will be zero;
 * that's why we start with p (kpivrw).
 
 * While we do this product, we also zero out the p'th row.
 */
static void c_ekketju_aux(COIN_REGISTER2 EKKfactinfo * COIN_RESTRICT2 fact,int sparse,
                  double * COIN_RESTRICT dluval, int * COIN_RESTRICT hrowi,
                  const int * COIN_RESTRICT mcstrt, const int * COIN_RESTRICT hpivco,
                  double * COIN_RESTRICT dwork1,
                  int *ipivp, int jpivrw, int stop_col)
{
  int ipiv = *ipivp;
  if (1&&ipiv<stop_col&&c_ekk_IsSet(fact->bitArray,ipiv)) {
    /* slack */
    int lastSlack = fact->lastSlack;
    int firstDo=hpivco[lastSlack];
    while (1) {
      assert (c_ekk_IsSet(fact->bitArray,ipiv));
      dwork1[ipiv] = -dwork1[ipiv];
      ipiv=hpivco[ipiv];      /* next column - generally ipiv+1 */
      if (ipiv==firstDo||ipiv>=stop_col)
      break;
    }
  }
  
  while(ipiv<stop_col) {
    double dv = dwork1[ipiv];
    int kx = mcstrt[ipiv];
    int nel = hrowi[kx];
    double dpiv = dluval[kx];
    int kcs = kx + 1;
    int kce = kx + nel;
    int iel;
    
    for (iel = kcs; iel <= kce; ++iel) {
      int irow = hrowi[iel];
      irow = UNSHIFT_INDEX(irow);
      dv -= dwork1[irow] * dluval[iel];
      if (irow == jpivrw) {
      break;
      }
    }
    
    /* assuming the p'th row is sparse,
     * this branch will be infrequently taken */
    if (iel <= kce) {
      int irow = hrowi[iel];
      /* irow == jpivrw */
      dv += dluval[iel];
      
      if (sparse) {
      /* delete this entry by overwriting it with the last */
      --nel;
      hrowi[kx] = nel;
      hrowi[iel] = hrowi[kce];
      dluval[iel] = dluval[kce];
      kce--;
      } else {
      /* we can't delete an entry from a dense column,
       * so we just zero it out */
      dluval[iel]=0.0;
      iel++;
      }
      
      /* finish up the remaining entries; same as above loop, but no check */
      for (; iel <= kce; ++iel) {
      irow = UNSHIFT_INDEX(hrowi[iel]);
      dv -= dwork1[irow] * dluval[iel];
      }
    }
    dwork1[ipiv] = dv * dpiv; /* divide by pivot */
    ipiv=hpivco[ipiv];  /* next column - generally ipiv+1 */
  }
  
  /* ? is it guaranteed that ipiv==stop_col at this point?? */
  *ipivp = ipiv;
}

/* dwork1 is assumed to be zeroed on entry */
static void c_ekketju(COIN_REGISTER EKKfactinfo * COIN_RESTRICT2 fact,double *dluval, int *hrowi,
                const int * COIN_RESTRICT mcstrt, const int * COIN_RESTRICT hpivco,
                double * COIN_RESTRICT dwork1,
                int kpivrw, int first_dense , int last_dense)
{
  int ipiv = hpivco[kpivrw];
  int jpivrw = SHIFT_INDEX(kpivrw);
  
  const int nrow  = fact->nrow;
  
  if (first_dense < last_dense &&
      mcstrt[ipiv] <= mcstrt[last_dense]) {
    /* There are dense columns, and
     * some dense columns precede the pivot column */
    
    /* first do any sparse columns "on the left" */
    c_ekketju_aux(fact, true, dluval, hrowi, mcstrt, hpivco, dwork1,
            &ipiv, jpivrw, first_dense);
    
    /* then do dense columns */
    c_ekketju_aux(fact, false, dluval, hrowi, mcstrt, hpivco, dwork1,
            &ipiv, jpivrw, last_dense+1);
    
    /* final sparse columns "on the right" ...*/
  } 
  /* ...are the same as sparse columns if there are no dense */
  c_ekketju_aux(fact, true, dluval, hrowi, mcstrt, hpivco, dwork1,
            &ipiv, jpivrw, nrow+1);
} /* c_ekketju */












/*#define PRINT_DEBUG*/
/* dwork1 is assumed to be zeroed on entry */
int c_ekketsj(COIN_REGISTER2 /*const*/ EKKfactinfo * COIN_RESTRICT2 fact,
          double * COIN_RESTRICT dwork1,
          int * COIN_RESTRICT mpt2, double dalpha, int orig_nincol,
          int npivot, int *nuspikp,
          const int ipivrw,int * spare)
{
  int nuspik      = *nuspikp;
  
  int * COIN_RESTRICT mpermu=fact->mpermu;
  
  int * COIN_RESTRICT hcoli   = fact->xecadr;
  double * COIN_RESTRICT dluval     = fact->xeeadr;
  int * COIN_RESTRICT mrstrt  = fact->xrsadr;
  int * COIN_RESTRICT hrowi   = fact->xeradr;
  int * COIN_RESTRICT mcstrt  = fact->xcsadr;
  int * COIN_RESTRICT hinrow  = fact->xrnadr;
  /*int *hincol   = fact->xcnadr;
    int *hpivro   = fact->krpadr;*/
  int * COIN_RESTRICT hpivco  = fact->kcpadr;
  double * COIN_RESTRICT de2val     = fact->xe2adr;
  
  const int nrow  = fact->nrow;
  const int ifRowCopy   = fact->rows_ok;
  
  int i, j=-1, k, i1, i2, k1;
  int kc, iel;
  double del3;
  int nroom;
  bool ifrows= (mrstrt[1] != 0);
  int kpivrw, jpivrw;
  int first_dense_mcstrt,last_dense_mcstrt;
  int nnentl;                 /* includes row stuff */
  int doSparse=(fact->if_sparse_update>0);
  
  int mcstrt_piv;
  int nincol=0;
  int * COIN_RESTRICT hpivco_new=fact->kcpadr+1;
  int * COIN_RESTRICT back=fact->back;
  int irtcod = 0;
  
  /* Parameter adjustments */
  de2val--;
  
  /* Function Body */
  if (!ifRowCopy) {
    doSparse=0;
    fact->if_sparse_update=-abs(fact->if_sparse_update);
  }
  if (npivot==1) {
    fact->num_resets=0;
  }
  kpivrw = mpermu[ipivrw];
#ifndef NDEBUG
  ets_count++;
  if (ets_check>=0&&ets_count>=ets_check) {
    printf("trouble\n");
  }
#endif
  mcstrt_piv=mcstrt[kpivrw];
  /* ndenuc - top has number deleted */
  if (fact->ndenuc) {
    first_dense_mcstrt = mcstrt[fact->first_dense];
    last_dense_mcstrt  = mcstrt[fact->last_dense];
  } else {
    first_dense_mcstrt=0;
    last_dense_mcstrt=0;
  }
  {
    int kdnspt = fact->nnetas - fact->nnentl;
    
    i1 = ((kdnspt - 1) + fact->R_etas_start[fact->nR_etas + 1]);
    /*i1 = -99999999;*/
    
    /* fact->R_etas_start[fact->nR_etas + 1] is -(the number of els in R) */
    nnentl = fact->nnetas - ((kdnspt - 1) + fact->R_etas_start[fact->nR_etas + 1]);
  }
  fact->demark=fact->nnentu+nnentl;
  jpivrw = SHIFT_INDEX(kpivrw);
  
  if (nuspik < 0) {
    goto L7000;
  } else if (nuspik == 0) {
    del3 = 0.;
  } else {
    del3 = 0.;
    i1 = fact->nnentu + 1;
    i2 = fact->nnentu + nuspik;
    if (fact->sortedEta) {
      /* binary search */
      if (hrowi[i2] == jpivrw) {
      /* sitting right on the end - easy */
      del3 = dluval[i2];
      --nuspik;
      } else {
      bool foundit = true;
      
      /* binary search - sort of implies hrowi is sorted */
      i = i1;
      if (hrowi[i] != jpivrw) {
        while (1) {
          i = (i1 + i2) >>1;
          if (i == i1) {
            foundit = false;
            break;
          }
          if (hrowi[i] < jpivrw) {
            i1 = i;
          } else if (hrowi[i] > jpivrw) {
            i2 = i;
          }
          else
            break;
        }
      }
      /* ??? what if we didn't find it? */
      
      if (foundit) {
        del3 = dluval[i];
        --nuspik;
        /* remove it and move the last element into its place */
        hrowi[i] = hrowi[nuspik + fact->nnentu+1];
        dluval[i] = dluval[nuspik + fact->nnentu+1];
      }
      }
    } else {
      /* search */
      for (i=i1;i<=i2;i++) {
      if (hrowi[i] == jpivrw) {
        del3 = dluval[i];
        --nuspik;
        /* remove it and move the last element into its place */
        hrowi[i] = hrowi[i2];
        dluval[i] = dluval[i2];
        break;
      }
      }
    }
  }
  
  /*      OLD COLUMN POINTERS */
  /* **************************************************************** */
  if (!ifRowCopy) {
    /*       old method */
    /*       DO U */
    c_ekketju(fact,dluval, hrowi, mcstrt, hpivco_new,
          dwork1, kpivrw,fact->first_dense,
          fact->last_dense);
  } else {
    
    /*       could take out of old column but lets try being crude */
    /*       try taking out */
    if (fact->xe2adr != 0&&doSparse) {
      
      /*
       * There is both a column and row representation of U.
       * For each row in the kpivrw'th column of the col U rep,
       * find its position in the U row rep and remove it
       * by overwriting it with the last element.
       */
      int k1x = mcstrt[kpivrw];
      int nel = hrowi[k1x];   /* yes, this is the nel, for the pivot */
      int k2x = k1x + nel;
      
      for (k = k1x + 1; k <= k2x; ++k) {
      int irow = UNSHIFT_INDEX(hrowi[k]);
      int kx = mrstrt[irow];
      int nel = hinrow[irow]-1;
      hinrow[irow]=nel;
      
      //int iel = c_ekkfrst2(&hcoli[kx], nel, kpivrw);
      //assert (iel>=0);
      //iel += kx;
      //int iel;
      int jlast = kx + nel;
      for (int iel=kx;iel<jlast;iel++) {
        if (kpivrw==hcoli[iel]) {
          hcoli[iel] = hcoli[jlast];
          de2val[iel] = de2val[jlast];
          break;
        }
      }
      //assert (iel<=jlast);
      }
    } else if (ifRowCopy) {
      /* still take out */
      int k1x = mcstrt[kpivrw];
      int nel = hrowi[k1x];   /* yes, this is the nel, for the pivot */
      int k2x = k1x + nel;
      
      for (k = k1x + 1; k <= k2x; ++k) {
      int irow = UNSHIFT_INDEX(hrowi[k]);
      int kx = mrstrt[irow];
      int nel = hinrow[irow]-1;
      hinrow[irow]=nel;
      int jlast = kx + nel ;
      for (;kx<jlast;kx++) {
        if (kpivrw==hcoli[kx]) {
          hcoli[kx] = hcoli[jlast];
          break;
        }
      }
      }
    }
    
    /*       add to row version */
    /* the updated column (alpha_p) was written to entries
     * nnentu+1..nnentu+nuspik by routine c_ekkftrn_ft.
     * That was just an intermediate value of the usual ftrn.
     */
    i1 = fact->nnentu + 1;
    i2 = fact->nnentu + nuspik;
    int * COIN_RESTRICT eta_last=mpermu+nrow*2+3;
    int * COIN_RESTRICT eta_next=eta_last+nrow+2;
    if (fact->xe2adr == 0||!doSparse) {
      /* we have column indices by row, but not the actual values */
      for (iel = i1; iel <= i2; ++iel) {
      int irow = UNSHIFT_INDEX(hrowi[iel]);
      int iput = hinrow[irow];
      int kput = mrstrt[irow];
      int nextRow=eta_next[irow];
      assert (kput>0);
      kput += iput;
      if (kput < mrstrt[nextRow]) {
        /* there is room - append the pivot column;
         * this corresponds making alpha_p the rightmost column of U (p. 268)*/
        hinrow[irow] = iput + 1;
        hcoli[kput] = kpivrw;
      } else {
        /* no room - switch off */
        doSparse=0;
        /* possible kpivrw 1 */
        k1 = mrstrt[kpivrw];
        mrstrt[1]=-1;
        fact->rows_ok = false;
        goto L1226;
      }
      }
    } else {
      if (! doSparse) {
      /* we have both column indices and values by row */
      /* just like loop above, but with extra assign to de2val */
      for (iel = i1; iel <= i2; ++iel) {
        int irow = UNSHIFT_INDEX(hrowi[iel]);
        int iput = hinrow[irow];
        int kput = mrstrt[irow];
        int nextRow=eta_next[irow];
        assert (kput>0);
        kput += iput;
        if (kput < mrstrt[nextRow]) {
          hinrow[irow] = iput + 1;
          hcoli[kput] = kpivrw;
          de2val[kput] = dluval[iel];
        } else {
          /* no room - switch off */
          doSparse=0;
          /* possible kpivrw 1 */
          k1 = mrstrt[kpivrw];
          mrstrt[1]=-1;
          fact->rows_ok = false;
          goto L1226;
        }
      }
      } else {
      for (iel = i1; iel <= i2; ++iel) {
        int j,k;
        int irow = UNSHIFT_INDEX(hrowi[iel]);
        int iput = hinrow[irow];
        k=mrstrt[irow]+iput;
        j=eta_next[irow];
        if (k >= mrstrt[j]) {
          /* no room - can we make some? */
          int klast=eta_last[nrow+1];
          int jput=mrstrt[klast]+hinrow[klast]+2;
          int distance=mrstrt[nrow+1]-jput;
          if (iput+1<distance) {
            /* this presumably copies the row to the end */
            int jn,jl;
            int kstart=mrstrt[irow];
            int nin=hinrow[irow];
            /* out */
            jn=eta_next[irow];
            jl=eta_last[irow];
            eta_next[jl]=jn;
            eta_last[jn]=jl;
            /* in */
            eta_next[klast]=irow;
            eta_last[nrow+1]=irow;
            eta_last[irow]=klast;
            eta_next[irow]=nrow+1;
            mrstrt[irow]=jput;
#if 0
            memcpy(&hcoli[jput],&hcoli[kstart],nin*sizeof(int));
            memcpy(&de2val[jput],&de2val[kstart],nin*sizeof(double));
#else
            c_ekkscpy(nin,hcoli+kstart,hcoli+jput);
            c_ekkdcpy(nin,
                  (de2val+kstart),(de2val+jput));
#endif
            k=jput+iput;
          } else {
            /* shuffle down */
            int spare=(fact->nnetas-fact->nnentu-fact->nnentl);
            if (spare>nrow<<1) {
            /* presumbly, this compacts the rows */
            int jrow,jput;
            if (1) {
              if (fact->num_resets<1000000) {
                int etasize =CoinMax(4*fact->nnentu+
                               (fact->nnetas-fact->nnentl)+1000,fact->eta_size);
                if (ifrows) {
                  fact->num_resets++;
                  if (npivot>40&&fact->num_resets<<4>npivot) {
                  fact->eta_size=static_cast<int>(1.05*fact->eta_size);
                  fact->num_resets=1000000;
                  }
                } else {
                  fact->eta_size=static_cast<int>(1.1*fact->eta_size);
                  fact->num_resets=1000000;
                }
                fact->eta_size=CoinMin(fact->eta_size,etasize);
                if (fact->maxNNetas>0&&fact->eta_size>
                  fact->maxNNetas) {
                  fact->eta_size=fact->maxNNetas;
                }
              }
            }
            jrow=eta_next[0];
            jput=1;
            for (j=0;j<nrow;j++) {
              int k,nin=hinrow[jrow];
              k=mrstrt[jrow];
              mrstrt[jrow]=jput;
              for (;nin;nin--) {
                hcoli[jput]=hcoli[k];
                de2val[jput++]=de2val[k++];
              }
              jrow=eta_next[jrow];
            }
            if (spare>nrow<<3) {
              spare=3;
            } else if (spare>nrow<<2) {
              spare=1;
            } else {
              spare=0;
            } 
            jput+=nrow*spare;;
            jrow=eta_last[nrow+1];
            for (j=0;j<nrow;j++) {
              int k,nin=hinrow[jrow];
              k=mrstrt[jrow]+nin;
              jput-=spare;
              for (;nin;nin--) {
                hcoli[--jput]=hcoli[--k];
                de2val[jput]=de2val[k];
              }
              mrstrt[jrow]=jput;
              jrow=eta_last[jrow];
            }
            /* set up for copy below */
            k=mrstrt[irow]+iput;
            } else {
            /* no room - switch off */
            doSparse=0;
            /* possible kpivrw 1 */
            k1 = mrstrt[kpivrw];
            mrstrt[1]=-1;
            fact->rows_ok = false;
            goto L1226;
            }
          }
        }
        /* now we have room - append the new value */
        hinrow[irow] = iput + 1;
        hcoli[k] = kpivrw;
        de2val[k] = dluval[iel];
      }
      }
    }
    
    /*       TAKE OUT ALL ELEMENTS IN PIVOT ROW */
    k1 = mrstrt[kpivrw];
    
  L1226:
    {
      int k2 = k1 + hinrow[kpivrw] - 1;
      
      /* "delete" the row */
      hinrow[kpivrw] = 0;
      j = 0;
      if (doSparse) {
      /* remove pivot row entries from the corresponding columns */
      for (k = k1; k <= k2; ++k) {
        int icol = hcoli[k];
        int kx = mcstrt[icol];
        int nel = hrowi[kx];
        for (iel = kx + 1; iel <= kx+nel; iel ++) {
          if (hrowi[iel] == jpivrw) {
            break;
          }
        }
        if (iel <= kx+nel) {
          /* this has to happen, right?? */
          
          /* copy the element into a temporary */
          dwork1[icol] = dluval[iel];
          mpt2[nincol++]=icol;
          /*nonzero[icol-1]=1;*/
          
          hrowi[kx]=nel-1;    /* column is shorter by one */
          j=1;
          hrowi[iel]=hrowi[kx+nel];
          dluval[iel]=dluval[kx+nel];
        }
      }
      if (j != 0) {
        /* now compute r', the new R transform */
        orig_nincol=c_ekkbtju_sparse(fact, dwork1,
                             mpt2, nincol,
                             spare);
        dwork1[kpivrw]=0.0;
      }
      } else {
      /* row version isn't ok (?) */
      for (k = k1; k <= k2; ++k) {
        int icol = hcoli[k];
        int kx = mcstrt[icol];
        int nel = hrowi[kx];
        j = kx+nel;
        int iel;
        for (iel=kx+1;iel<=j;iel++) {
          if (hrowi[iel]==jpivrw)
            break;
        }
        dwork1[icol] = dluval[iel];
        if (kx<first_dense_mcstrt || kx>last_dense_mcstrt) {
          hrowi[kx] = nel - 1;      /* shorten column */
          /* not packing - presumably column isn't sorted */
          hrowi[iel] = hrowi[j];
          dluval[iel] = dluval[j];
        } else {
          /* dense element - just zero it */
          dluval[iel]=0.0;
        }
      }
      if (j != 0) {
        /* Find first nonzero */
        int ipiv = hpivco_new[kpivrw];
        while(ipiv<=nrow) {
          if (!dwork1[ipiv]) {
            ipiv=hpivco_new[ipiv];
          } else {
            break;
          }
        }
        if (ipiv<=nrow) {
          /*       DO U */
          /* now compute r', the new R transform */
          c_ekkbtju(fact, dwork1,
                  ipiv);
        }
      }
      }
    }
  }
  
  if (kpivrw==fact->first_dense) {
    /* increase until valid pivot */
    fact->first_dense=hpivco_new[fact->first_dense];
  } else if (kpivrw==fact->last_dense) {
    fact->last_dense=back[fact->last_dense];
  }
  if (fact->first_dense==fact->last_dense) {
    fact->ndenuc=0;
    fact->first_dense=0;
    fact->last_dense=-1;
  }
  
  if (! (ifRowCopy && j==0)) {
    
    /*     increase amount of work on Etas */
    
    /* **************************************************************** */
    /*       DO ROW ETAS IN L */
    {
      if (!doSparse) {
      dwork1[kpivrw] = 0.;
#if 0 
      orig_nincol=c_ekksczr(fact,nrow, dwork1, mpt2);
      del3=c_ekkputl(fact, mpt2, dwork1, del3, 
                 orig_nincol, nuspik);
#else
      orig_nincol=c_ekkputl2(fact,
                  dwork1, &del3, 
                 nuspik);
#endif
      } else {
      del3=c_ekkputl(fact, mpt2,
                  dwork1, del3, 
                 orig_nincol, nuspik);
      }
    }
    if (orig_nincol != 0) {
      /* STORE AS A ROW VECTOR */
      int n = fact->nR_etas+1;
      int i1 = fact->R_etas_start[n];
      fact->nR_etas=n;
      fact->R_etas_start[n + 1] = i1 - orig_nincol;
      hpivco[fact->nR_etas + nrow+3] = kpivrw;
    }
  }
  
  /*       CHECK DEL3 AGAINST DALPHA/DOUT */
  {
    int kx = mcstrt[kpivrw];
    double dout = dluval[kx];
    double dcheck = fabs(dalpha / dout);
    double difference=0.0;
    if (fabs(del3) > CoinMin(1.0e-8,fact->drtpiv*0.99999)) {
      double checkTolerance;
      if ( fact->npivots < 2 ) {
      checkTolerance = 1.0e-5;
      } else if ( fact->npivots < 10 ) {
      checkTolerance = 1.0e-6;
      } else if ( fact->npivots < 50 ) {
      checkTolerance = 1.0e-8;
      } else {
      checkTolerance = 1.0e-9;
      }
      difference = fabs(1.0-fabs(del3)/dcheck);
      if (difference > 0.1*checkTolerance) {
      if (difference < checkTolerance||
          (difference<1.0e-7&&fact->npivots>=50)) {
        irtcod=1;
#ifdef PRINT_DEBUG
        printf("mildly bad %g after %d pivots, etsj %g ftncheck %g ftnalpha %g\n",
             difference,fact->npivots,del3,dcheck,dalpha);
#endif    
      } else {
        irtcod=2;
#ifdef PRINT_DEBUG
        printf("bad %g after %d pivots, etsj %g ftncheck %g ftnalpha %g\n",
             difference,fact->npivots,del3,dcheck,dalpha);
#endif    
      }
      }
    } else {
      irtcod=2;
#ifdef PRINT_DEBUG
      printf("bad small %g after %d pivots, etsj %g ftncheck %g ftnalpha %g\n",
           difference,fact->npivots,del3,dcheck,dalpha);
#endif    
    }
    if (irtcod>1)
      goto L8000;
    fact->npivots++;
  }
  
  mcstrt[kpivrw] = fact->nnentu;
  dluval[fact->nnentu] = 1. / del3; /* new pivot */
  hrowi[fact->nnentu] = nuspik;     /* new nelems */
#ifndef NDEBUG
  {
    int lastSlack = fact->lastSlack;
    int firstDo=hpivco_new[lastSlack];
    int ipiv=hpivco_new[0];
    int now = fact->numberSlacks;
    if (now) {
      while (1) {
      if (ipiv>fact->nrow||ipiv==firstDo)
        break;
      assert (c_ekk_IsSet(fact->bitArray,ipiv));
      ipiv=hpivco_new[ipiv];
      }
      if (ipiv<=fact->nrow) {
      while (1) {
        if (ipiv>fact->nrow)
          break;
        assert (!c_ekk_IsSet(fact->bitArray,ipiv));
        ipiv=hpivco_new[ipiv];
      }
      }
    }
  }
#endif
  {
    /* do new hpivco */
    int inext=hpivco_new[kpivrw];
    int iback=back[kpivrw];
    if (inext!=nrow+1) {
      int ilast=back[nrow+1];
      hpivco_new[iback]=inext;
      back[inext]=iback;
      assert (hpivco_new[ilast]==nrow+1);
      hpivco_new[ilast]=kpivrw;
      back[kpivrw]=ilast;
      hpivco_new[kpivrw]=nrow+1;
      back[nrow+1]=kpivrw;
    }
  }
  {
    int lastSlack = fact->lastSlack;
    int now = fact->numberSlacks;
    if (now&&mcstrt_piv<=mcstrt[lastSlack]) {
      if (c_ekk_IsSet(fact->bitArray,kpivrw)) {
      /*printf("piv %d lastSlack %d\n",mcstrt_piv,lastSlack);*/
      fact->numberSlacks--;   
      now--;
      /* one less slack */
      c_ekk_Unset(fact->bitArray,kpivrw);
      if (now&&kpivrw==lastSlack) {
        int i;
        int ipiv,jpiv;
        ipiv=hpivco_new[0];
        for (i=0;i<now-1;i++)
          ipiv=hpivco_new[ipiv];
        lastSlack=ipiv;
        jpiv=hpivco_new[ipiv];
        assert (c_ekk_IsSet(fact->bitArray,ipiv));
        assert (!c_ekk_IsSet(fact->bitArray,jpiv)||jpiv>fact->nrow);
        fact->lastSlack = lastSlack;
      } else if (!now) {
        fact->lastSlack=0;
      }
      }
    }
    fact->firstNonSlack=hpivco_new[lastSlack];
#ifndef NDEBUG
    {
      int lastSlack = fact->lastSlack;
      int firstDo=hpivco_new[lastSlack];
      int ipiv=hpivco_new[0];
      int now = fact->numberSlacks;
      if (now) {
      while (1) {
        if (ipiv>fact->nrow||ipiv==firstDo)
          break;
        assert (c_ekk_IsSet(fact->bitArray,ipiv));
        ipiv=hpivco_new[ipiv];
      }
      if (ipiv<=fact->nrow) {
        while (1) {
          if (ipiv>fact->nrow)
            break;
          assert (!c_ekk_IsSet(fact->bitArray,ipiv));
          ipiv=hpivco_new[ipiv];
        }
      }
      }
    }
#endif
  }
  
  fact->nnentu += nuspik;
  {
    int kdnspt = fact->nnetas - fact->nnentl;
    
    /* fact->R_etas_start[fact->nR_etas + 1] is -(the number of els in R) */
    nnentl = fact->nnetas - ((kdnspt - 1) + fact->R_etas_start[fact->nR_etas + 1]);
  }
  fact->demark = (fact->nnentu + nnentl) - fact->demark;
  
  /*     if need to redo row version */
  if (! fact->rows_ok&&fact->first_dense>=fact->last_dense) {
    int extraSpace=10000;
    int spareSpace;
    if (fact->if_sparse_update>0) {
      spareSpace=(fact->nnetas-fact->nnentu-fact->nnentl);
    } else {
      /* missing out nnentl stuff */
      spareSpace=fact->nnetas-fact->nnentu;
    }
    /*       save clean row copy if enough room */
    nroom = spareSpace / nrow;
    
    if ((fact->nnentu<<3)>150*fact->maxinv) {
      extraSpace=150*fact->maxinv;
    } else {
      extraSpace=fact->nnentu<<3;
    }
    
    ifrows = false;
    if (fact->nnetas>fact->nnentu+fact->nnentl+extraSpace) {
      ifrows = true;
    }
    if (nroom < 5) {
      ifrows = false;
    }
    
    if (nroom > CoinMin(50, fact->maxinv - (fact->iterno - fact->iterin))) {
      ifrows = true;
    }
    
#ifdef PRINT_DEBUGx
    printf(" redoing row copy %d %d %d\n",ifrows,nroom,spareSpace);
#endif
    if (1) {
      if (fact->num_resets<1000000) {
      if (ifrows) {
        fact->num_resets++;
        if (npivot>40&&fact->num_resets<<4>npivot) {
          fact->eta_size=static_cast<int>(1.05*fact->eta_size);
          fact->num_resets=1000000;
        }
      } else {
        fact->eta_size=static_cast<int>(1.1*fact->eta_size);
        fact->num_resets=1000000;
      }
      if (fact->maxNNetas>0&&fact->eta_size>
          fact->maxNNetas) {
        fact->eta_size=fact->maxNNetas;
      }
      }
    }
    fact->rows_ok = ifrows;
    if (ifrows) {
      int ibase = 1;
      c_ekkizero(nrow,&hinrow[1]);
      for (i = 1; i <= nrow; ++i) {
      int kx = mcstrt[i];
      int nel = hrowi[kx];
      int kcs = kx + 1;
      int kce = kx + nel;
      for (kc = kcs; kc <= kce; ++kc) {
        int irow = UNSHIFT_INDEX(hrowi[kc]);
        if (dluval[kc]) {
          hinrow[irow]++;
        }
      }
      }
      int * eta_last=mpermu+nrow*2+3;
      int * eta_next=eta_last+nrow+2;
      eta_next[0]=1;
      for (i = 1; i <= nrow; ++i) {
      eta_next[i]=i+1;
      eta_last[i]=i-1;
      mrstrt[i] = ibase;
      ibase = ibase + hinrow[i] + nroom;
      hinrow[i] = 0;
      }
      eta_last[nrow+1]=nrow;
      //eta_next[nrow+1]=nrow+2;
      mrstrt[nrow+1]=ibase;
      if (fact->xe2adr == 0) {
      for (i = 1; i <= nrow; ++i) {
        int kx = mcstrt[i];
        int nel = hrowi[kx];
        int kcs = kx + 1;
        int kce = kx + nel;
        for (kc = kcs; kc <= kce; ++kc) {
          if (dluval[kc]) {
            int irow = UNSHIFT_INDEX(hrowi[kc]);
            int iput = hinrow[irow];
            assert (irow);
            hcoli[mrstrt[irow] + iput] = i;
            hinrow[irow] = iput + 1;
          }
        }
      }
      } else {
      for (i = 1; i <= nrow; ++i) {
        int kx = mcstrt[i];
        int nel = hrowi[kx];
        int kcs = kx + 1;
        int kce = kx + nel;
        for (kc = kcs; kc <= kce; ++kc) {
          int irow = UNSHIFT_INDEX(hrowi[kc]);
          int iput = hinrow[irow];
          hcoli[mrstrt[irow] + iput] = i;
          de2val[mrstrt[irow] + iput] = dluval[kc];
          hinrow[irow] = iput + 1;
        }
      }
      }
    } else {
      mrstrt[1] = 0;
      if (fact->if_sparse_update>0&&fact->iterno-fact->iterin>100) {
      goto L7000;
      }
    }
  }
  goto L8000;
  
  /*       OUT OF SPACE - COULD PACK DOWN */
 L7000:
  irtcod = 1;
#ifdef PRINT_DEBUG
  printf(" out of space\n");
#endif    
  if (1) {
    if ((npivot<<3)<fact->nbfinv) {
      /* low on space */
      if (npivot<10) {
      fact->eta_size=fact->eta_size<<1;
      } else {
      double ratio=fact->nbfinv;
      double ratio2=npivot<<3;
      ratio=ratio/ratio2;
      if (ratio>2.0) {
        ratio=2.0;
      } /* endif */
      fact->eta_size=static_cast<int>(ratio*fact->eta_size);
      } /* endif */
    } else {
      fact->eta_size=static_cast<int>(1.05*fact->eta_size);
    } /* endif */
    if (fact->maxNNetas>0&&fact->eta_size>
      fact->maxNNetas) {
      fact->eta_size=fact->maxNNetas;
    }
  }
  
  /* ================= IF ERROR SHOULD WE GET RID OF LAST ITERATION??? */
 L8000:
  
  *nuspikp = nuspik;
  return (irtcod);
} /* c_ekketsj */
static void c_ekkftj4p(COIN_REGISTER2 const EKKfactinfo * COIN_RESTRICT2 fact, 
                 double * COIN_RESTRICT dwork1, int firstNonZero)
{
  /* this is where the L factors start, because this is the place
   * where c_ekktria starts laying them down (see initialization of xnetal).
   */
  int lstart=fact->lstart;
  const int * COIN_RESTRICT hpivco  = fact->kcpadr;
  int firstLRow = hpivco[lstart];
  if (firstNonZero>firstLRow) {
    lstart += firstNonZero-firstLRow;
  }
  assert (firstLRow==fact->firstLRow);
  int jpiv=hpivco[lstart];
  const double * COIN_RESTRICT dluval     = fact->xeeadr;
  const int * COIN_RESTRICT hrowi   = fact->xeradr;
  const int * COIN_RESTRICT mcstrt  = fact->xcsadr+lstart;
  int ndo=fact->xnetal-lstart;
  int i, iel;
  
  /* find first non-zero */
  for (i=0;i<ndo;i++) {
    if (dwork1[i+jpiv]!=0.0)
      break;
  }
  for (; i < ndo; ++i) {
    double dv = dwork1[i+jpiv];
    
    if (dv != 0.) {
      int kce1 = mcstrt[i + 1] ;
      
      for (iel = mcstrt[i]; iel > kce1; --iel) {
      int irow0 = hrowi[iel];
      SHIFT_REF(dwork1, irow0) += dv * dluval[iel];
      }
    }
  }
  
} /* c_ekkftj4p */

/*
 * This version is more efficient for input columns that are sparse.
 * It is instructive to consider the case of an especially sparse column,
 * which is a slack.  The slack for row r has exactly one non-zero element,
 * in row r, which is +-1.0.  Let pr = mpermu[r].
 * In this case, nincol==1 and mpt[0] == pr on entry.
 * if mpt[0] == pr <= jpiv
 * then this slack is completely unaffected by L;
 *    this is reflected by the fact that save_where = last
 *    after the first loop, so none of the remaining loops
 *    ever execute,
 * else if mpt[0] == pr > jpiv, but pr-jpiv > ndo
 * then the slack is also unaffected by L, this time because
 *    its row is "after" L.  During factorization, it may
 *    be the case that the first part of the basis is upper
 *    triangular (c_ekktria), but it may also be the case that the
 *    last part of the basis is upper triangular (in which case the
 *    L triangle gets "chopped off" on the right).  In both cases,
 *    no L entries are required.  Since in this case the tests
 *    (i<=ndo) will fail (and dwork1[ipiv]==1.0), the code will
 *    do nothing.
 * else if mpt[0] == pr > jpiv and pr-jpiv <= ndo
 * then the slack *is* affected by L.
 *    the for-loop inside the second while-loop will discover
 *    that none of the factors for the corresponding column of L
 *    are non-zero in the slack column, so last will not be incremented.
 *    We multiply the eta-vector, and the last loop does nothing.
 */
static int c_ekkftj4_sparse(COIN_REGISTER2 const EKKfactinfo * COIN_RESTRICT2 fact,
               double * COIN_RESTRICT dwork1, int * COIN_RESTRICT mpt,
               int nincol,int * COIN_RESTRICT spare)
{
  const int nrow  = fact->nrow;
  /* this is where the L factors start, because this is the place
   * where c_ekktria starts laying them down (see initialization of xnetal).
   */
  int lstart=fact->lstart;
  const int * COIN_RESTRICT hpivco  = fact->kcpadr;
  const double * COIN_RESTRICT dluval     = fact->xeeadr;
  const int * COIN_RESTRICT hrowi   = fact->xeradr;
  const int * COIN_RESTRICT mcstrt  = fact->xcsadr+lstart-1;
  double tolerance = fact->zeroTolerance;
  int jpiv=hpivco[lstart]-1;
  char * COIN_RESTRICT nonzero=fact->nonzero;
  int ndo=fact->xnetalval;
  int k,nStack;
  int nList=0;
  int iPivot;
  int * COIN_RESTRICT list = spare;
  int * COIN_RESTRICT stack = spare+nrow;
  int * COIN_RESTRICT next = stack+nrow;
  double dv;
  int iel;
  int nput=0,kput=nrow;
  int check=jpiv+ndo+1;
  const int * COIN_RESTRICT mcstrt2 = mcstrt-jpiv;
  
#if 0
  for (k=0;k<nincol;k++) {
    iPivot=mpt[k];
    nonzero[iPivot]=0;
  }
#endif
  for (k=0;k<nincol;k++) {
    nStack=1;
    iPivot=mpt[k];
    if (nonzero[iPivot]!=1&&iPivot>jpiv&&iPivot<check) {
      stack[0]=iPivot;
      next[0]=mcstrt2[iPivot+1]+1;
      while (nStack) {
      int kPivot,j;
      /* take off stack */
      kPivot=stack[--nStack];
      if (nonzero[kPivot]!=1&&kPivot>jpiv&&kPivot<check) {
        j=next[nStack];
        if (j>mcstrt2[kPivot]) {
          /* finished so mark */
          list[nList++]=kPivot;
          nonzero[kPivot]=1;
        } else {
          kPivot=UNSHIFT_INDEX(hrowi[j]);
          /* put back on stack */
          next[nStack++] ++;
          if (!nonzero[kPivot]) {
            /* and new one */
            stack[nStack]=kPivot;
            nonzero[kPivot]=2;
            next[nStack++]=mcstrt2[kPivot+1]+1;
          }
        }
      } else if (kPivot>=check) {
        list[--kput]=kPivot;
        nonzero[kPivot]=1;
      }
      }
    } else if (nonzero[iPivot]!=1) {
      /* nothing to do (except check size at end) */
      list[--kput]=iPivot;
      nonzero[iPivot]=1;
    }
  }
  for (k=nList-1;k>=0;k--) {
    double dv;
    iPivot = list[k];
    dv = dwork1[iPivot];
    nonzero[iPivot]=0;
    if (fabs(dv) > tolerance) {
      /* the same code as in c_ekkftj4p */
      int kce1 = mcstrt2[iPivot + 1];
      for (iel = mcstrt2[iPivot]; iel > kce1; --iel) {
      int irow0 = hrowi[iel];
      SHIFT_REF(dwork1, irow0) += dv * dluval[iel];
      }
      mpt[nput++]=iPivot;
    } else {
      dwork1[iPivot]=0.0;     /* force to zero, not just near zero */
    }
  }
  /* check remainder */
  for (k=kput;k<nrow;k++) {
    iPivot = list[k];
    nonzero[iPivot]=0;
    dv = dwork1[iPivot];
    if (fabs(dv) > tolerance) {
      mpt[nput++]=iPivot;
    } else {
      dwork1[iPivot]=0.0;     /* force to zero, not just near zero */
    }
  }
  
  return (nput);
} /* c_ekkftj4 */
/*
 * This applies the R transformations of the F-T LU update procedure,
 * equation 3.11 on p. 270 in the 1972 Math Programming paper.
 * Note that since the non-zero off-diagonal elements are in a row,
 * multiplying an R by a column is a reduction, not like applying
 * L or U.
 *
 * Note that this may introduce new non-zeros in dwork1,
 * since an hpivco entry may correspond to a zero element,
 * and that some non-zeros in dwork1 may be cancelled.
 */
static int c_ekkftjl_sparse3(COIN_REGISTER2 const EKKfactinfo * COIN_RESTRICT2 fact,
                double * COIN_RESTRICT dwork1, 
                int * COIN_RESTRICT mpt,
                int * COIN_RESTRICT hput, double * COIN_RESTRICT dluput ,
                int nincol)
{
  int i;
  int knext;
  int ipiv;
  double dv;
  const double * COIN_RESTRICT dluval = fact->R_etas_element+1;
  const int * COIN_RESTRICT hrowi = fact->R_etas_index+1;
  const int * COIN_RESTRICT mcstrt = fact->R_etas_start;
  int ndo=fact->nR_etas;
  double tolerance = fact->zeroTolerance;
  const int * COIN_RESTRICT hpivco = fact->hpivcoR;
  /* and make cleaner */
  hput++;
  dluput++;
  
  /* DO ANY ROW TRANSFORMATIONS */
  
  /* Function Body */
  /* mpt has correct list of nonzeros */
  if (ndo != 0) {
    knext = mcstrt[1];
    for (i = 1; i <= ndo; ++i) {
      int k1 = knext;   /* == mcstrt[i] */
      int iel;
      ipiv = hpivco[i];
      dv = dwork1[ipiv];
      /* put on list if not there */
      if (!dv)
      mpt[nincol++]=ipiv;
      knext = mcstrt[i + 1];
      
      for (iel = knext ; iel < k1; ++iel) {
      int irow = hrowi[iel];
      dv += SHIFT_REF(dwork1, irow) * dluval[iel];
      }
      /* (1) if dwork[ipiv] == 0.0, then this may add a non-zero.
       * (2) if dwork[ipiv] != 0.0, then this may cancel out a non-zero.
       */
      if (fabs(dv) > tolerance) {
      dwork1[ipiv]=dv;
      } else {
      dwork1[ipiv] = 1.0e-128;
      }
    }
  }
  knext=0;
  for (i=0; i<nincol; i++) {
    ipiv=mpt[i];
    dv=dwork1[ipiv];
    if (fabs(dv) > tolerance) {
      hput[knext]=SHIFT_INDEX(ipiv);
      dluput[knext]=dv;
      mpt[knext++]=ipiv;
    } else {
      dwork1[ipiv]=0.0;
    }
  }
  return knext;
} /* c_ekkftjl */

static int c_ekkftjl_sparse2(COIN_REGISTER2 const EKKfactinfo * COIN_RESTRICT2 fact,
                double * COIN_RESTRICT dwork1, 
                int * COIN_RESTRICT mpt,
                int nincol)
{
  double tolerance = fact->zeroTolerance;
  const double * COIN_RESTRICT dluval = fact->R_etas_element+1;
  const int * COIN_RESTRICT hrowi = fact->R_etas_index+1;
  const int * COIN_RESTRICT mcstrt = fact->R_etas_start;
  int ndo=fact->nR_etas;
  const int * COIN_RESTRICT hpivco = fact->hpivcoR;
  int i;
  int knext;
  int ipiv;
  double dv;
  
  /* DO ANY ROW TRANSFORMATIONS */
  
  /* Function Body */
  /* mpt has correct list of nonzeros */
  if (ndo != 0) {
    knext = mcstrt[1];
    for (i = 1; i <= ndo; ++i) {
      int k1 = knext;   /* == mcstrt[i] */
      int iel;
      ipiv = hpivco[i];
      dv = dwork1[ipiv];
      /* put on list if not there */
      if (!dv)
      mpt[nincol++]=ipiv;
      knext = mcstrt[i + 1];
      
      for (iel = knext ; iel < k1; ++iel) {
      int irow = hrowi[iel];
      dv += SHIFT_REF(dwork1, irow) * dluval[iel];
      }
      /* (1) if dwork[ipiv] == 0.0, then this may add a non-zero.
       * (2) if dwork[ipiv] != 0.0, then this may cancel out a non-zero.
       */
      if (fabs(dv) > tolerance) {
      dwork1[ipiv]=dv;
      } else {
      dwork1[ipiv] = 1.0e-128;
      }
    }
  }
  knext=0;
  for (i=0; i<nincol; i++) {
    ipiv=mpt[i];
    dv=dwork1[ipiv];
    if (fabs(dv) > tolerance) {
      mpt[knext++]=ipiv;
    } else {
      dwork1[ipiv]=0.0;
    }
  }
  return knext;
} /* c_ekkftjl */

static void c_ekkftjl(COIN_REGISTER2 const EKKfactinfo * COIN_RESTRICT2 fact,
           double * COIN_RESTRICT dwork1)
{
  double tolerance = fact->zeroTolerance;
  const double * COIN_RESTRICT dluval = fact->R_etas_element+1;
  const int * COIN_RESTRICT hrowi = fact->R_etas_index+1;
  const int * COIN_RESTRICT mcstrt = fact->R_etas_start;
  int ndo=fact->nR_etas;
  const int * COIN_RESTRICT hpivco = fact->hpivcoR;
  int i;
  int knext;
  
  /* DO ANY ROW TRANSFORMATIONS */
  
  /* Function Body */
  if (ndo != 0) {
    /*
     * The following three lines are here just to ensure that this
     * new formulation of the loop has exactly the same effect
     * as the original.
     */
    {
      int ipiv = hpivco[1];
      double dv = dwork1[ipiv];
      dwork1[ipiv] = (fabs(dv) > tolerance) ? dv : 0.0;
    }
    
    knext = mcstrt[1];
    for (i = 1; i <= ndo; ++i) {
      int k1 = knext;   /* == mcstrt[i] */
      int ipiv = hpivco[i];
      double dv = dwork1[ipiv];
      int iel;
      //#define UNROLL3 2
#ifndef UNROLL3
#if CLP_OSL==2||CLP_OSL==3
#define UNROLL3 2
#else
#define UNROLL3 1
#endif
#endif
      knext = mcstrt[i + 1];
      
#if UNROLL3<2
      for (iel = knext ; iel < k1; ++iel) {
      int irow = hrowi[iel];
      dv += SHIFT_REF(dwork1, irow) * dluval[iel];
      }
#else
      iel = knext;
      if (((k1-knext)&1)!=0) {
      int irow = hrowi[iel];
      dv += SHIFT_REF(dwork1, irow) * dluval[iel];
      iel++;
      }
      for ( ; iel < k1; iel+=2) {
      int irow0 = hrowi[iel];
      double dval0 = dluval[iel];
      int irow1 = hrowi[iel+1];
      double dval1 = dluval[iel+1];
      dv += SHIFT_REF(dwork1, irow0) * dval0;
      dv += SHIFT_REF(dwork1, irow1) * dval1;
      }
#endif
      /* (1) if dwork[ipiv] == 0.0, then this may add a non-zero.
       * (2) if dwork[ipiv] != 0.0, then this may cancel out a non-zero.
       */
      dwork1[ipiv] = (fabs(dv) > tolerance) ? dv : 0.0;
    }
  }
} /* c_ekkftjl */
/* this assumes it is ok to reference back[loop_limit] */
/* another 3 seconds from a ~570 second run can be trimmed
 * by using two routines, one with scan==true and the other false,
 * since that eliminates the branch instructions involving them
 * entirely.  This was how the code was originally written.
 * However, I'm still hoping that eventually we can use
 * C++ templates to do that for us automatically.
 */
static void
c_ekkftjup_scan_aux(COIN_REGISTER2 const EKKfactinfo * COIN_RESTRICT2 fact,
              double * COIN_RESTRICT dwork1, double * COIN_RESTRICT dworko ,
              int loop_limit, int *ip, int ** mptp)
{
  const double * COIN_RESTRICT dluval     = fact->xeeadr+1;
  const int * COIN_RESTRICT hrowi   = fact->xeradr+1;
  const int * COIN_RESTRICT mcstrt  = fact->xcsadr;
  const int * COIN_RESTRICT hpivro  = fact->krpadr;
  const int * COIN_RESTRICT back=fact->back;
  double tolerance = fact->zeroTolerance;
  int ipiv = *ip;
  double dv = dwork1[ipiv];
  
  int * mptX = *mptp;
  assert (mptX);
  while (ipiv != loop_limit) {
    int next_ipiv = back[ipiv];
    
    dwork1[ipiv] = 0.0;
#ifndef UNROLL4
#define UNROLL4 2
#endif
    /* invariant:  dv == dwork1[ipiv] */
    
    /* in the case of world.mps with dual, this condition is true
     * only 20-60% of the time. */
    if (fabs(dv) > tolerance) {
      const int kx = mcstrt[ipiv];
      const int nel = hrowi[kx-1];
      const double dpiv = dluval[kx-1];
#if UNROLL4>1
      const int * hrowi2=hrowi+kx;
      const int * hrowi2end=hrowi2+nel;
      const double * dluval2=dluval+kx;
#else
      int iel;
#endif
      
      dv*=dpiv;
      
      /*
       * The following loop is the unrolled version of this:
       *
       * for (iel = kx+1; iel <= kx + nel; iel++) {
       *   SHIFT_REF(dwork1, hrowi[iel]) -= dv * dluval[iel];
       * }
       */
#if UNROLL4<2
      iel = kx;
      if (nel&1) {
      int irow = hrowi[iel];
      double dval=dluval[iel];
      SHIFT_REF(dwork1, irow) -= dv*dval;
      iel++;
      }
      for (; iel < kx + nel; iel+=2) {
      int irow0 = hrowi[iel];
      int irow1 = hrowi[iel+1];
      double dval0=dluval[iel];
      double dval1=dluval[iel+1];
      double d0=SHIFT_REF(dwork1, irow0);
      double d1=SHIFT_REF(dwork1, irow1);
      
      d0-=dv*dval0;
      d1-=dv*dval1;
      SHIFT_REF(dwork1, irow0)=d0;
      SHIFT_REF(dwork1, irow1)=d1;
      } /* end loop */
#else
      if ((nel&1)!=0) {
      int irow = *hrowi2;
      double dval=*dluval2;
      SHIFT_REF(dwork1, irow) -= dv*dval;
      hrowi2++;
      dluval2++;
      }
      for (; hrowi2 < hrowi2end; hrowi2 +=2,dluval2 +=2) {
      int irow0 = hrowi2[0];
      int irow1 = hrowi2[1];
      double dval0=dluval2[0];
      double dval1=dluval2[1];
      double d0=SHIFT_REF(dwork1, irow0);
      double d1=SHIFT_REF(dwork1, irow1);
      
      d0-=dv*dval0;
      d1-=dv*dval1;
      SHIFT_REF(dwork1, irow0)=d0;
      SHIFT_REF(dwork1, irow1)=d1;
      }
#endif
      /* put this down here so that dv is less likely to cause a stall */
      if (fabs(dv) >= tolerance) {
      int iput=hpivro[ipiv];
      dworko[iput]=dv;
      *mptX++=iput-1;
      }
    }
    
    dv = dwork1[next_ipiv];
    ipiv=next_ipiv;
  } /* endwhile */
  
  *mptp = mptX;
  *ip = ipiv;
}
static void c_ekkftjup_aux3(COIN_REGISTER2 const EKKfactinfo * COIN_RESTRICT2 fact,
                       double * COIN_RESTRICT dwork1, double * COIN_RESTRICT dworko,
                    const int * COIN_RESTRICT back,
                    const int * COIN_RESTRICT hpivro,
                    int *ipivp, int loop_limit,
                    int **mptXp)
  
{
  double tolerance = fact->zeroTolerance;
  int ipiv = *ipivp;
  if (ipiv!=loop_limit) {
    int *mptX = *mptXp;
    
    double dv = dwork1[ipiv];
    
    do {
      int next_ipiv = back[ipiv];
      double next_dv=dwork1[next_ipiv];
      
      dwork1[ipiv]=0.0;
      
      if (fabs(dv)>=tolerance) {
      int iput=hpivro[ipiv];
      dworko[iput]=dv;
      *mptX++=iput-1;
      }
      
      ipiv = next_ipiv;
      dv = next_dv;
    } while (ipiv!=loop_limit);
    
    *mptXp = mptX;
    *ipivp = ipiv;
  }
}
static void c_ekkftju_dense(const double *dluval, 
                      const int * COIN_RESTRICT hrowi, 
                      const int * COIN_RESTRICT mcstrt,
                      const int * COIN_RESTRICT back,
                      double * COIN_RESTRICT dwork1, 
                      int * start, int last,
               int offset , double *densew)
{
  int ipiv=*start;
  
  while (ipiv>last ) {
    const int ipiv1=ipiv;
    double dv1=dwork1[ipiv1];
    ipiv=back[ipiv];
    if (fabs(dv1) > 1.0e-14) {
      const int kx1 = mcstrt[ipiv1];
      const int nel1 = hrowi[kx1-1];
      const double dpiv1 = dluval[kx1-1];
      
      int iel,k;
      const int n1=offset+ipiv1;    /* number in dense part */
      
      const int nsparse1=nel1-n1;
      const int k1=kx1+nsparse1;
      const double *dlu1=&dluval[k1];
      
      int ipiv2=back[ipiv1];
      const int nskip=ipiv1-ipiv2;
      
      dv1*=dpiv1;
      
      dwork1[ipiv1]=dv1;
      
      for (k = n1 - (nskip-1) -1; k >=0 ; k--) {
        const double dval = dv1*dlu1[k];
      double dv2=densew[k]-dval;
        ipiv=back[ipiv];
        if (fabs(dv2) > 1.0e-14) {
          const int kx2 = mcstrt[ipiv2];
          const int nel2 = hrowi[kx2-1];
          const double dpiv2 = dluval[kx2-1];
        
          /* number in dense part is k */
          const int nsparse2=nel2-k;
        
          const int k2=kx2+nsparse2;
          const double *dlu2=&dluval[k2];
        
          dv2*=dpiv2;
          densew[k]=dv2;      /* was dwork1[ipiv2]=dv2; */
        
          k--;
        
        /*
         * The following loop is the unrolled version of:
         *
         * for (; k >= 0; k--) {
         *   densew[k]-=dv1*dlu1[k]+dv2*dlu2[k];
         * }
         */
          if ((k&1)==0) {
            densew[k]-=dv1*dlu1[k]+dv2*dlu2[k];
            k--;
          }
          for (; k >=0 ; k-=2) {
            double da,db;
            da=densew[k];
            db=densew[k-1];
            da-=dv1*dlu1[k];
            db-=dv1*dlu1[k-1];
            da-=dv2*dlu2[k];
            db-=dv2*dlu2[k-1];
            densew[k]=da;
            densew[k-1]=db;
          }
        /* end loop */
        
        /*
         * The following loop is the unrolled version of:
         *
         * for (iel=kx2+nsparse2-1; iel >= kx2; iel--) {
         *   SHIFT_REF(dwork1, hrowi[iel]) -= dv2*dluval[iel];
         * }
         */
        iel=kx2+nsparse2-1;
          if ((nsparse2&1)!=0) {
            int irow0 = hrowi[iel];
            double dval=dluval[iel];
            SHIFT_REF(dwork1,irow0) -= dv2*dval;
            iel--;
          }
          for (; iel >=kx2 ; iel-=2) {
            double dval0 = dluval[iel];
          double dval1 = dluval[iel-1];
            int irow0 = hrowi[iel];
            int irow1 = hrowi[iel-1];
            double d0 = SHIFT_REF(dwork1, irow0);
            double d1 = SHIFT_REF(dwork1, irow1);
          
            d0-=dv2*dval0;
            d1-=dv2*dval1;
          SHIFT_REF(dwork1, irow0) = d0;
          SHIFT_REF(dwork1, irow1) = d1;
          }
        /* end loop */
        
        } else {
          densew[k]=0.0;
          /* skip if next deleted */
          k-=ipiv2-ipiv-1;
          ipiv2=ipiv;
          if (ipiv<last) {
          k--;
          for (; k >=0 ; k--) {
            double dval;
            dval=dv1*dlu1[k];
            densew[k]=densew[k]-dval;
          }
        }
        }
      }
      
      /*
       * The following loop is the unrolled version of:
       *
       * for (iel=kx1+nsparse1-1; iel >= kx1; iel--) {
       *   SHIFT_REF(dwork1, hrowi[iel]) -= dv1*dluval[iel];
       * }
       */
      iel=kx1+nsparse1-1;
      if ((nsparse1&1)!=0) {
        int irow0 = hrowi[iel];
        double dval=dluval[iel];
        SHIFT_REF(dwork1, irow0) -= dv1*dval;
        iel--;
      }
      for (; iel >=kx1 ; iel-=2) {
        double dval0=dluval[iel];
        double dval1=dluval[iel-1];
        int irow0 = hrowi[iel];
        int irow1 = hrowi[iel-1];
        double d0=SHIFT_REF(dwork1, irow0);
        double d1=SHIFT_REF(dwork1, irow1);
      
        d0-=dv1*dval0;
        d1-=dv1*dval1;
        SHIFT_REF(dwork1, irow0) = d0;
        SHIFT_REF(dwork1, irow1) = d1;
      }
      /* end loop */
    } else {
      dwork1[ipiv1]=0.0;
    } /* endif */
  } /* endwhile */
  *start=ipiv;
}

/* do not use return value if mpt==0 */
/* using dual, this is usually called via c_ekkftrn_ft, from c_ekksdul
 * (so mpt is non-null).
 * it is generally called every iteration, but sometimes several iterations
 * are skipped (null moves?).
 *
 * generally, back[i] == i-1 (initialized in c_ekkshfv towards the end).
 */
static int c_ekkftjup(COIN_REGISTER3 const EKKfactinfo * COIN_RESTRICT2 fact,
           double * COIN_RESTRICT dwork1, int last,
           double * COIN_RESTRICT dworko , int * COIN_RESTRICT mpt)
{
  const double * COIN_RESTRICT dluval     = fact->xeeadr;
  const int * COIN_RESTRICT hrowi   = fact->xeradr;
  const int * COIN_RESTRICT mcstrt  = fact->xcsadr;
  const int * COIN_RESTRICT hpivro  = fact->krpadr;
  double tolerance = fact->zeroTolerance;
  int ndenuc=fact->ndenuc;
  const int first_dense=fact->first_dense;
  const int last_dense=fact->last_dense;
  int i;
  int * mptX = mpt;
  
  const int nrow        = fact->nrow;
  const int * COIN_RESTRICT back=fact->back;
  int ipiv=back[nrow+1];
  
  if (last_dense>first_dense&&mcstrt[ipiv]>=mcstrt[last_dense]) {
    c_ekkftjup_scan_aux(fact,
                  dwork1, dworko, last_dense, &ipiv,
                  &mptX);
    
    {
      int j;
      int n=0;
      const int firstDense    = nrow- ndenuc+1;
      double *densew = &dwork1[firstDense];
      int offset;
      
      /* check first dense to see where in triangle it is */
      int last=first_dense;
      const int k1=mcstrt[last];
      const int k2=k1+hrowi[k1];
      
      for (j=k2; j>k1; j--) {
        int irow = UNSHIFT_INDEX(hrowi[j]);
        if (irow<firstDense) {
          break;
        } else {
#ifdef DEBUG
          if (irow!=last-1) {
            abort();
          }
#endif
          last=irow;
          n++;
        }
      }
      offset=n-first_dense;
      i=ipiv;
      /* loop counter i may be modified by this call */
      c_ekkftju_dense(&dluval[1],&hrowi[1],mcstrt,back,
                dwork1, &i, first_dense,offset,densew);
      
      c_ekkftjup_aux3(fact,dwork1, dworko, back, hpivro, &ipiv, i, &mptX);
    }
  }
  
  c_ekkftjup_scan_aux(fact,
                dwork1, dworko, last, &ipiv,
                &mptX);
  
  if (ipiv!=0) {
    double dv = dwork1[ipiv];
    
    do {
      int next_ipiv = back[ipiv];
      double next_dv=dwork1[next_ipiv];
      
      dwork1[ipiv]=0.0;
      
      if (fabs(dv)>=tolerance) {
      int iput=hpivro[ipiv];
      dworko[iput]=-dv;
      *mptX++=iput-1;
      }
      
      ipiv = next_ipiv;
      dv = next_dv;
    } while (ipiv!=0);
    
  }
  return static_cast<int>(mptX-mpt);
}
/* this assumes it is ok to reference back[loop_limit] */
/* another 3 seconds from a ~570 second run can be trimmed
 * by using two routines, one with scan==true and the other false,
 * since that eliminates the branch instructions involving them
 * entirely.  This was how the code was originally written.
 * However, I'm still hoping that eventually we can use
 * C++ templates to do that for us automatically.
 */
static void
c_ekkftjup_scan_aux_pack(COIN_REGISTER2 const EKKfactinfo * COIN_RESTRICT2 fact,
                   double * COIN_RESTRICT dwork1, double * COIN_RESTRICT dworko ,
                   int loop_limit, int *ip, int ** mptp)
{
  double tolerance = fact->zeroTolerance;
  const double *dluval  = fact->xeeadr+1;
  const int *hrowi      = fact->xeradr+1;
  const int *mcstrt     = fact->xcsadr;
  const int *hpivro     = fact->krpadr;
  const int * back=fact->back;
  int ipiv = *ip;
  double dv = dwork1[ipiv];
  
  int * mptX = *mptp;
#if 0
  int inSlacks=0;
  int lastSlack;
  if (fact->numberSlacks!=0)
    lastSlack=fact->lastSlack;
  else
    lastSlack=0;
  if (c_ekk_IsSet(fact->bitArray,ipiv)) {
    printf("already in slacks - ipiv %d\n",ipiv);
    inSlacks=1;
    return;
  }
#endif
  assert (mptX);
  while (ipiv != loop_limit) {
    int next_ipiv = back[ipiv];
#if 0
    if (ipiv==lastSlack) {
      printf("now in slacks - ipiv %d\n",ipiv);
      inSlacks=1;
      break;
    }
    if (inSlacks) {
      assert (c_ekk_IsSet(fact->bitArray,ipiv));
      assert (dluval[mcstrt[ipiv]-1]==-1.0);
      assert (hrowi[mcstrt[ipiv]-1]==0);
    }
#endif
    dwork1[ipiv] = 0.0;
    /* invariant:  dv == dwork1[ipiv] */
    
    /* in the case of world.mps with dual, this condition is true
     * only 20-60% of the time. */
    if (fabs(dv) > tolerance) {
      const int kx = mcstrt[ipiv];
      const int nel = hrowi[kx-1];
      const double dpiv = dluval[kx-1];
#ifndef UNROLL5
#define UNROLL5 2
#endif
#if UNROLL5>1
      const int * hrowi2=hrowi+kx;
      const int * hrowi2end=hrowi2+nel;
      const double * dluval2=dluval+kx;
#else
      int iel;
#endif
      
      dv*=dpiv;
      
      /*
       * The following loop is the unrolled version of this:
       *
       * for (iel = kx+1; iel <= kx + nel; iel++) {
       *   SHIFT_REF(dwork1, hrowi[iel]) -= dv * dluval[iel];
       * }
       */
#if UNROLL5<2
      iel = kx;
      if (nel&1) {
      int irow = hrowi[iel];
      double dval=dluval[iel];
      SHIFT_REF(dwork1, irow) -= dv*dval;
      iel++;
      }
      for (; iel < kx + nel; iel+=2) {
      int irow0 = hrowi[iel];
      int irow1 = hrowi[iel+1];
      double dval0=dluval[iel];
      double dval1=dluval[iel+1];
      double d0=SHIFT_REF(dwork1, irow0);
      double d1=SHIFT_REF(dwork1, irow1);
      
      d0-=dv*dval0;
      d1-=dv*dval1;
      SHIFT_REF(dwork1, irow0)=d0;
      SHIFT_REF(dwork1, irow1)=d1;
      } /* end loop */
#else
      if ((nel&1)!=0) {
      int irow = *hrowi2;
      double dval=*dluval2;
      SHIFT_REF(dwork1, irow) -= dv*dval;
      hrowi2++;
      dluval2++;
      }
      for (; hrowi2 < hrowi2end; hrowi2 +=2,dluval2 +=2) {
      int irow0 = hrowi2[0];
      int irow1 = hrowi2[1];
      double dval0=dluval2[0];
      double dval1=dluval2[1];
      double d0=SHIFT_REF(dwork1, irow0);
      double d1=SHIFT_REF(dwork1, irow1);
      
      d0-=dv*dval0;
      d1-=dv*dval1;
      SHIFT_REF(dwork1, irow0)=d0;
      SHIFT_REF(dwork1, irow1)=d1;
      }
#endif
      /* put this down here so that dv is less likely to cause a stall */
      if (fabs(dv) >= tolerance) {
      int iput=hpivro[ipiv];
      *dworko++=dv;
      *mptX++=iput-1;
      }
    }
    
    dv = dwork1[next_ipiv];
    ipiv=next_ipiv;
  } /* endwhile */
  
  *mptp = mptX;
  *ip = ipiv;
}
static void c_ekkftjup_aux3_pack(COIN_REGISTER2 const EKKfactinfo * COIN_RESTRICT2 fact,
                          double * COIN_RESTRICT dwork1, double * COIN_RESTRICT dworko,
                         const int * COIN_RESTRICT back,
                         const int * COIN_RESTRICT hpivro,
                         int *ipivp, int loop_limit,
                         int **mptXp)
  
{
  double tolerance = fact->zeroTolerance;
  int ipiv = *ipivp;
  if (ipiv!=loop_limit) {
    int *mptX = *mptXp;
    
    double dv = dwork1[ipiv];
    do {
      int next_ipiv = back[ipiv];
      double next_dv=dwork1[next_ipiv];
      
      dwork1[ipiv]=0.0;
      
      if (fabs(dv)>=tolerance) {
      int iput=hpivro[ipiv];
      *dworko++=dv;
      *mptX++=iput-1;
      }
      
      ipiv = next_ipiv;
      dv = next_dv;
    } while (ipiv!=loop_limit);
    
    *mptXp = mptX;
    
    *ipivp = ipiv;
  }
}

/* do not use return value if mpt==0 */
/* using dual, this is usually called via c_ekkftrn_ft, from c_ekksdul
 * (so mpt is non-null).
 * it is generally called every iteration, but sometimes several iterations
 * are skipped (null moves?).
 *
 * generally, back[i] == i-1 (initialized in c_ekkshfv towards the end).
 */
static int c_ekkftjup_pack(COIN_REGISTER3 const EKKfactinfo * COIN_RESTRICT2 fact,
              double * COIN_RESTRICT dwork1, int last,
              double * COIN_RESTRICT dworko , int * COIN_RESTRICT mpt)
{
  const double * COIN_RESTRICT dluval     = fact->xeeadr;
  const int * COIN_RESTRICT hrowi   = fact->xeradr;
  const int * COIN_RESTRICT mcstrt  = fact->xcsadr;
  const int * COIN_RESTRICT hpivro  = fact->krpadr;
  double tolerance = fact->zeroTolerance;
  int ndenuc=fact->ndenuc;
  const int first_dense=fact->first_dense;
  const int last_dense=fact->last_dense;
  int * mptX = mpt; int * mptY = mpt;
  
  const int nrow        = fact->nrow;
  const int * COIN_RESTRICT back=fact->back;
  int ipiv=back[nrow+1];
  assert (mpt);
  
  if (last_dense>first_dense&&mcstrt[ipiv]>=mcstrt[last_dense]) {
    c_ekkftjup_scan_aux_pack(fact,
                     dwork1, dworko, last_dense, &ipiv,
                     &mptX );
    /* adjust */
    dworko+= (mptX-mpt);
    mpt=mptX;
    {
      int j;
      int n=0;
      const int firstDense    = nrow- ndenuc+1;
      double *densew = &dwork1[firstDense];
      int offset;
      
      /* check first dense to see where in triangle it is */
      int last=first_dense;
      const int k1=mcstrt[last];
      const int k2=k1+hrowi[k1];
      
      for (j=k2; j>k1; j--) {
        int irow = UNSHIFT_INDEX(hrowi[j]);
        if (irow<firstDense) {
          break;
        } else {
#ifdef DEBUG
          if (irow!=last-1) {
            abort();
          }
#endif
          last=irow;
          n++;
        }
      }
      offset=n-first_dense;
      int ipiv2=ipiv;
      /* loop counter i may be modified by this call */
      c_ekkftju_dense(&dluval[1],&hrowi[1],mcstrt,back,
                dwork1, &ipiv2, first_dense,offset,densew);
      
      c_ekkftjup_aux3_pack(fact,dwork1, dworko, back, hpivro, &ipiv, ipiv2,&mptX);
      /* adjust dworko */
      dworko += (mptX-mpt);
      mpt=mptX;
    }
  }
  
  c_ekkftjup_scan_aux_pack(fact,
                   dwork1, dworko, last, &ipiv,
                   &mptX );
  /* adjust dworko */
  dworko += (mptX-mpt);
  while (ipiv!=0) {
    double dv = dwork1[ipiv];
    int next_ipiv = back[ipiv];
    
    dwork1[ipiv]=0.0;
    
    if (fabs(dv)>=tolerance) {
      int iput=hpivro[ipiv];
      *dworko++=-dv;
      *mptX++=iput-1;
    }
    
    ipiv = next_ipiv;
  }
  
  return static_cast<int>(mptX-mptY);
}
static int c_ekkftju_sparse_a(COIN_REGISTER2 const EKKfactinfo * COIN_RESTRICT2 fact,
                      int * COIN_RESTRICT mpt,
                      int nincol,int * COIN_RESTRICT spare)
{
  const int * COIN_RESTRICT hrowi   = fact->xeradr+1;
  const int * COIN_RESTRICT mcstrt  = fact->xcsadr;
  const int nrow        = fact->nrow;
  char * COIN_RESTRICT nonzero=fact->nonzero;
  
  int k,nStack,kx,nel;
  int nList=0;
  int iPivot;
  /*int kkk=nincol;*/
  int * COIN_RESTRICT list = spare;
  int * COIN_RESTRICT stack = spare+nrow;
  int * COIN_RESTRICT next = stack+nrow;
  for (k=0;k<nincol;k++) {
    nStack=1;
    iPivot=mpt[k];
    stack[0]=iPivot;
    next[0]=0;
    while (nStack) {
      int kPivot,j;
      /* take off stack */
      kPivot=stack[--nStack];
      if (nonzero[kPivot]!=1) {
      kx = mcstrt[kPivot];
      nel = hrowi[kx-1];
      j=next[nStack];
      if (j==nel) {
        /* finished so mark */
        list[nList++]=kPivot;
        nonzero[kPivot]=1;
      } else {
        kPivot=hrowi[kx+j];
        /* put back on stack */
        next[nStack++] ++;
        if (!nonzero[kPivot]) {
          /* and new one */
          stack[nStack]=kPivot;
          nonzero[kPivot]=2;
          next[nStack++]=0;
          /*kkk++;*/
        }
      }
      }
    }
  }
  return (nList);
}
static int c_ekkftju_sparse_b(COIN_REGISTER2 const EKKfactinfo * COIN_RESTRICT2 fact,
                      double * COIN_RESTRICT dwork1, 
                      double * COIN_RESTRICT dworko , int * COIN_RESTRICT mpt,
                      int nList,int * COIN_RESTRICT spare)
{
  
  const double * COIN_RESTRICT dluval     = fact->xeeadr+1;
  const int * COIN_RESTRICT hrowi   = fact->xeradr+1;
  const int * COIN_RESTRICT mcstrt  = fact->xcsadr;
  const int * COIN_RESTRICT hpivro  = fact->krpadr;
  double tolerance = fact->zeroTolerance;
  char * COIN_RESTRICT nonzero=fact->nonzero;
  int i,k,kx,nel;
  int iPivot;
  /*int kkk=nincol;*/
  int * COIN_RESTRICT list = spare;
  i=nList-1;
  nList=0;
  for (;i>=0;i--) {
    double dpiv;
    double dv;
    iPivot = list[i];
    /*printf("pivot %d %d\n",i,iPivot);*/
    dv=dwork1[iPivot];
    kx = mcstrt[iPivot];
    nel = hrowi[kx-1];
    dwork1[iPivot]=0.0;
    dpiv = dluval[kx-1];
    dv*=dpiv;
    nonzero[iPivot]=0;
    iPivot=hpivro[iPivot];
    if (fabs(dv)>=tolerance) {
      *dworko++=dv;
      mpt[nList++]=iPivot-1;
      for (k = kx; k < kx+nel; k++) {
        double dval;
        double dd;
        int irow = hrowi[k];
        dval=dluval[k];
        dd=dwork1[irow];
        dd-=dv*dval;
        dwork1[irow]=dd;
      }
    }
  }
  return (nList);
}
/* dwork1 = (B^-1)dwork1;
 * I think dpermu[1..nrow+1] is zeroed on exit (?)
 * I don't think it is expected to have any particular value on entry (?)
 */
int c_ekkftrn(COIN_REGISTER const EKKfactinfo * COIN_RESTRICT2 fact, 
          double * COIN_RESTRICT dwork1, 
          double * COIN_RESTRICT dpermu, int * COIN_RESTRICT mpt,int numberNonZero)
{
  const int * COIN_RESTRICT mpermu = fact->mpermu;
  int lastNonZero;
  int firstNonZero = c_ekkshfpi_list2(mpermu+1, dwork1+1, dpermu, mpt,
                              numberNonZero,&lastNonZero);
  if (fact->nnentl&&lastNonZero>=fact->firstLRow) {
    /* dpermu = (L^-1)dpermu */
    c_ekkftj4p(fact, dpermu, firstNonZero);
  }
  
  
  int lastSlack;
  
  /* dpermu = (R^-1) dpermu */
  c_ekkftjl(fact, dpermu);
  
  assert (fact->numberSlacks!=0||!fact->lastSlack);
  lastSlack=fact->lastSlack;
  
  /* dwork1 = (U^-1)dpermu; dpermu zeroed (?) */
  return c_ekkftjup(fact, 
              dpermu, lastSlack, dwork1, mpt);
  
} /* c_ekkftrn */

int c_ekkftrn_ft(COIN_REGISTER EKKfactinfo * COIN_RESTRICT2 fact,
             double * COIN_RESTRICT dwork1_ft, int * COIN_RESTRICT mpt_ft, int *nincolp_ft)
{
  double * COIN_RESTRICT dpermu_ft = fact->kadrpm;
  int * COIN_RESTRICT spare = reinterpret_cast<int *>(fact->kp1adr);
  int nincol      = *nincolp_ft;
  
  int nuspik;
  double * COIN_RESTRICT dluvalPut = fact->xeeadr+fact->nnentu+1;
  int * COIN_RESTRICT hrowiPut      = fact->xeradr+fact->nnentu+1;
  
  const int nrow  = fact->nrow;
  /* mpermu contains the permutation */
  const int * COIN_RESTRICT mpermu=fact->mpermu;
  
  int lastSlack;
  
  int kdnspt = fact->nnetas - fact->nnentl;
  bool isRoom = (fact->nnentu + (nrow << 1) < (kdnspt - 1) 
             + fact->R_etas_start[fact->nR_etas + 1]);
  
  /* say F-T will be sorted */
  fact->sortedEta=1;
  
  assert (fact->numberSlacks!=0||!fact->lastSlack);
  lastSlack=fact->lastSlack;
  
  if (fact->if_sparse_update>0 &&  DENSE_THRESHOLD<nrow) {
    /* iterating so c_ekkgtcl will have list */
    {
      /* in order for this to make sense, nonzero[1..nrow] must already be zeroed */
      c_ekkshfpi_list3(mpermu+1, dwork1_ft, dpermu_ft, mpt_ft, nincol);
      
      /* it may be the case that the basis was entirely upper-triangular */
      if (fact->nnentl) {
      nincol = 
        c_ekkftj4_sparse(fact,
                   dpermu_ft,  mpt_ft,
                   nincol,spare);
      }
      /*       DO ROW ETAS IN L */
      if (isRoom) {
      ++fact->nnentu;
      nincol=
        c_ekkftjl_sparse3(fact,
                    dpermu_ft, 
                    mpt_ft, hrowiPut,
                    dluvalPut,nincol);
      nuspik = nincol;
        /* temporary */
      /* say not sorted */
      fact->sortedEta=0;
      } else {
      /* no room */
      nuspik=-3;
      nincol=
        c_ekkftjl_sparse2(fact,
                    dpermu_ft, 
                    mpt_ft,  nincol);
      }
      /*         DO U */
      if (DENSE_THRESHOLD>nrow-fact->numberSlacks) {
      nincol = c_ekkftjup_pack(fact, 
                         dpermu_ft,lastSlack, dwork1_ft, 
                         mpt_ft);
      } else {
      nincol= c_ekkftju_sparse_a(fact,
                         mpt_ft,
                         nincol, spare);
      nincol = c_ekkftju_sparse_b(fact,
                          dpermu_ft, 
                          dwork1_ft , mpt_ft,
                          nincol, spare);
      }
    }
  } else {
    int lastNonZero;
    int firstNonZero = c_ekkshfpi_list(mpermu+1, dwork1_ft, dpermu_ft, 
                               mpt_ft, nincol,&lastNonZero);
    if (fact->nnentl&&lastNonZero>=fact->firstLRow) {
      /* dpermu_ft = (L^-1)dpermu_ft */
      c_ekkftj4p(fact, dpermu_ft, firstNonZero);
    }
    
    
    /* dpermu_ft = (R^-1) dpermu_ft */
    c_ekkftjl(fact, dpermu_ft);
    
    if (isRoom) {
      
      /*        fake start to allow room for pivot */
      /* dluval[fact->nnentu...] = non-zeros of dpermu_ft;
       * hrowi[fact->nnentu..] = indices of these non-zeros;
       * near-zeros in dluval flattened
       */
      ++fact->nnentu;
      nincol= c_ekkscmv(fact,fact->nrow, dpermu_ft, hrowiPut,
                  dluvalPut);
      
      /*
       * note that this is not the value of nincol determined by c_ekkftjup. 
       * For Forrest-Tomlin update we want vector before U
       * this vector will replace one in U
       */
      nuspik = nincol;
    } else {
      /* no room */
      nuspik = -3;
    }
    
    /* dwork1_ft = (U^-1)dpermu_ft; dpermu_ft zeroed (?) */
    nincol = c_ekkftjup_pack(fact, 
                     dpermu_ft, lastSlack, dwork1_ft, mpt_ft);
  }
  
  *nincolp_ft = nincol;
  return (nuspik);
} /* c_ekkftrn */
void c_ekkftrn2(COIN_REGISTER EKKfactinfo * COIN_RESTRICT2 fact, double * COIN_RESTRICT dwork1,
           double * COIN_RESTRICT dpermu1,int * COIN_RESTRICT mpt1, int *nincolp,
           double * COIN_RESTRICT dwork1_ft, int * COIN_RESTRICT mpt_ft, int *nincolp_ft)
{
  double * COIN_RESTRICT dluvalPut = fact->xeeadr+fact->nnentu+1;
  int * COIN_RESTRICT hrowiPut      = fact->xeradr+fact->nnentu+1;
  
  const int nrow  = fact->nrow;
  /* mpermu contains the permutation */
  const int * COIN_RESTRICT mpermu=fact->mpermu;
  
  int lastSlack;
  assert (fact->numberSlacks!=0||!fact->lastSlack);
  lastSlack=fact->lastSlack;
  
  int nincol      = *nincolp_ft;

  /* using dwork1 instead double *dpermu_ft = fact->kadrpm; */
  int * spare = reinterpret_cast<int *>(fact->kp1adr);
  
  int kdnspt = fact->nnetas - fact->nnentl;
  bool isRoom = (fact->nnentu + (nrow << 1) < (kdnspt - 1) 
             + fact->R_etas_start[fact->nR_etas + 1]);
  
  /* say F-T will be sorted */
  fact->sortedEta=1;
  int lastNonZero;
  int firstNonZero = c_ekkshfpi_list2(mpermu+1, dwork1+1, dpermu1, 
                              mpt1, *nincolp,&lastNonZero);
  if (fact->nnentl&&lastNonZero>=fact->firstLRow) {
    /* dpermu1 = (L^-1)dpermu1 */
    c_ekkftj4p(fact, dpermu1, firstNonZero);
  }
  
  
  if (fact->if_sparse_update>0 &&  DENSE_THRESHOLD<nrow) {
    /* iterating so c_ekkgtcl will have list */
    {
      /* in order for this to make sense, nonzero[1..nrow] must already be zeroed */
      c_ekkshfpi_list3(mpermu+1, dwork1_ft, dwork1, mpt_ft, nincol);
      
      /* it may be the case that the basis was entirely upper-triangular */
      if (fact->nnentl) {
      nincol = 
        c_ekkftj4_sparse(fact,
                   dwork1, mpt_ft,
                   nincol,spare);
      }
      /*       DO ROW ETAS IN L */
      if (isRoom) {
      ++fact->nnentu;
      nincol=
        c_ekkftjl_sparse3(fact,
                    dwork1, 
                    mpt_ft, hrowiPut,
                    dluvalPut,
                    nincol);
      fact->nuspike = nincol;
      /* say not sorted */
      fact->sortedEta=0;
      } else {
      /* no room */
      fact->nuspike=-3;
      nincol=
        c_ekkftjl_sparse2(fact,
                    dwork1, 
                    mpt_ft, nincol);
      }
    }
  } else {
    int lastNonZero;
    int firstNonZero = c_ekkshfpi_list(mpermu+1, dwork1_ft, dwork1, 
                               mpt_ft, nincol,&lastNonZero);
    if (fact->nnentl&&lastNonZero>=fact->firstLRow) {
      /* dpermu_ft = (L^-1)dpermu_ft */
      c_ekkftj4p(fact, dwork1, firstNonZero);
    }
    
    
    /* dpermu_ft = (R^-1) dpermu_ft */
    c_ekkftjl(fact, dwork1);
    
    if (isRoom) {
      
      /*        fake start to allow room for pivot */
      /* dluval[fact->nnentu...] = non-zeros of dpermu_ft;
       * hrowi[fact->nnentu..] = indices of these non-zeros;
       * near-zeros in dluval flattened
       */
      ++fact->nnentu;
      nincol= c_ekkscmv(fact,fact->nrow, dwork1, 
                  hrowiPut,
                  dluvalPut);
      
      /*
       * note that this is not the value of nincol determined by c_ekkftjup. 
       * For Forrest-Tomlin update we want vector before U
       * this vector will replace one in U
       */
      fact->nuspike = nincol;
    } else {
      /* no room */
      fact->nuspike = -3;
    }
  }

  
  /* dpermu1 = (R^-1) dpermu1 */
  c_ekkftjl(fact, dpermu1);
  
  /*         DO U */
  if (fact->if_sparse_update<=0 ||  DENSE_THRESHOLD>nrow-fact->numberSlacks) {
    nincol = c_ekkftjup_pack(fact,
                     dwork1,lastSlack, dwork1_ft, mpt_ft);
  } else {
    nincol= c_ekkftju_sparse_a(fact,
                       mpt_ft,
                       nincol, spare);
    nincol = c_ekkftju_sparse_b(fact,
                        dwork1, 
                        dwork1_ft , mpt_ft,
                        nincol, spare);
  }
  *nincolp_ft = nincol;
  /* dwork1 = (U^-1)dpermu1; dpermu1 zeroed (?) */
  *nincolp = c_ekkftjup(fact,
                  dpermu1,lastSlack, dwork1, mpt1);

}



Generated by  Doxygen 1.6.0   Back to index