/* $Id: CoinOslFactorization2.cpp 1585 2013-04-06 20:42:02Z stefan $ */
/*
  Copyright (C) 1987, 2009, International Business Machines
  Corporation and others.  All Rights Reserved.

  This code is licensed under the terms of the Eclipse Public License (EPL).
*/
/*
  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 "CoinOslFactorization.hpp"
#include "CoinOslC.h"
#include "CoinFinite.hpp"

#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];
#ifdef CLP_REUSE_ETAS
	double temp=dluval[iel];
	dluval[iel] = dluval[kce];
	hrowi[kce]=jpivrw;
	dluval[kce]=temp;
#else
	dluval[iel] = dluval[kce];
#endif
	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);
#ifdef MORE_DEBUG
  {
    const int * COIN_RESTRICT hrowi = fact->R_etas_index;
    const int * COIN_RESTRICT mcstrt = fact->R_etas_start;
    int ndo=fact->nR_etas;
    int knext;
  
    knext = mcstrt[ndo + 1];
    for (int i = ndo; i > 0; --i) {
      int k1 = knext;
      knext = mcstrt[i];
      for (int j = k1+1; j < knext; j++) {
	assert (hrowi[j]>0&&hrowi[j]<100000);
      }
    }
  }
#endif
  
  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];
#if 0 //ndef 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);
  
#ifdef CLP_REUSE_ETAS
  double del3Orig=0.0;
#endif
  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;
	}
      }
    }
  }
#ifdef CLP_REUSE_ETAS
  del3Orig=del3;
#endif
  
  /*      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 jlast = kx + nel;
	for (int iel=kx;iel<jlast;iel++) {
	  if (kpivrw==hcoli[iel]) {
	    hcoli[iel] = hcoli[jlast];
	    de2val[iel] = de2val[jlast];
	    break;
	  }
	}
      }
    } 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-3);
	      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];
#ifdef CLP_REUSE_ETAS
	    hrowi[kx+nel]=jpivrw;
	    dluval[kx+nel]=dwork1[icol];
#endif
	  }
	}
	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];
#ifdef CLP_REUSE_ETAS
	    hrowi[j]=jpivrw;
	    dluval[j]=dwork1[icol];
#endif
	  } 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;
#ifdef CLP_REUSE_ETAS
  {
    int * putSeq = fact->xrsadr+2*fact->nrowmx+2;
    int * position = putSeq+fact->maxinv;
    int * putStart = position+fact->maxinv;
    putStart[fact->nrow+fact->npivots-1]=fact->nnentu;
  }
#endif
  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;
	  ipiv=hpivco_new[0];
	  for (i=0;i<now-1;i++)
	    ipiv=hpivco_new[ipiv];
	  lastSlack=ipiv;
	  assert (c_ekk_IsSet(fact->bitArray,ipiv));
	  assert (!c_ekk_IsSet(fact->bitArray,hpivco_new[ipiv])||hpivco_new[ipiv]>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;
#ifdef CLP_REUSE_ETAS
  if (fact->first_dense>=fact->last_dense) {
    // save
    fact->nnentu++;
    dluval[fact->nnentu]=del3Orig;
    hrowi[fact->nnentu]=kpivrw;
    int * putSeq = fact->xrsadr+2*fact->nrowmx+2;
    int * position = putSeq+fact->maxinv;
    int * putStart = position+fact->maxinv;
    int nnentu_at_factor=putStart[fact->nrow]&0x7fffffff;
    //putStart[fact->nrow+fact->npivots]=fact->nnentu+1;
    int where;
    if (mcstrt_piv<nnentu_at_factor) {
      // original LU
      where=kpivrw-1;
    } else {
      // could do binary search
      int * look = putStart+fact->nrow;
      for (where=fact->npivots-1;where>=0;where--) {
	if (mcstrt_piv==(look[where]&0x7fffffff))
	  break;
      }
      assert (where>=0);
      where += fact->nrow;
    }
    position[fact->npivots-1]=where;
    if (orig_nincol == 0) {
      // flag
      putStart[fact->nrow+fact->npivots-1] |= 0x80000000;
    }
  }
#endif
  {
    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;
#ifdef MORE_DEBUG
  for (int i=1;i<=fact->nrow;i++) {
    int kx=mcstrt[i];
    int nel=hrowi[kx];
    for (int j=0;j<nel;j++) {
      assert (i!=hrowi[j+kx+1]);
    }
  }
#endif
#ifdef CLP_REUSE_ETAS
  fact->save_nnentu=fact->nnentu;
#endif
  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;
  
  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];
      bool onList = (dv!=0.0);
      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 (onList) {
	if (fabs(dv) > tolerance) {
	  dwork1[ipiv]=dv;
	} else {
	  dwork1[ipiv] = 1.0e-128;
	}
      } else {
	if (fabs(dv) > tolerance) {
	  /* put on list if not there */
	  mpt[nincol++]=ipiv;
	  dwork1[ipiv]=dv;
	} 
      }
    }
  }
  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];
      bool onList = (dv!=0.0);
      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 (onList) {
	if (fabs(dv) > tolerance) {
	  dwork1[ipiv]=dv;
	} else {
	  dwork1[ipiv] = 1.0e-128;
	}
      } else {
	if (fabs(dv) > tolerance) {
	  /* put on list if not there */
	  mpt[nincol++]=ipiv;
	  dwork1[ipiv]=dv;
	} 
      }
    }
  }
  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 - 2) 
		 + 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;
#ifdef CLP_REUSE_ETAS
  bool skipStuff = (fact->reintro>=0);
  
  int save_nR_etas=fact->nR_etas;
  int * save_hpivcoR=fact->hpivcoR;
  int * save_R_etas_start=fact->R_etas_start;
  if (skipStuff) {
    // just move
    int * putSeq = fact->xrsadr+2*fact->nrowmx+2;
    int * position = putSeq+fact->maxinv;
    int * putStart = position+fact->maxinv;
    memset(dwork1_ft,0,nincol*sizeof(double));
    int iPiv=fact->reintro;
    int start=putStart[iPiv]&0x7fffffff;
    int end=putStart[iPiv+1]&0x7fffffff;
    double * COIN_RESTRICT dluval	= fact->xeeadr;
    int * COIN_RESTRICT hrowi	= fact->xeradr;
    double dValue;
    if (fact->reintro<fact->nrow) {
      iPiv++;
      dValue=1.0/dluval[start++];
    } else {
      iPiv=hrowi[--end];
      dValue=dluval[end];
      start++;
      int ndoSkip=0;
      for (int i=fact->nrow;i<fact->reintro;i++) {
	if ((putStart[i]&0x80000000)==0)
	ndoSkip++;
      }
      fact->nR_etas-=ndoSkip;
      fact->hpivcoR+=ndoSkip;
      fact->R_etas_start+=ndoSkip;
    }
    dpermu_ft[iPiv]=dValue;
    if (fact->if_sparse_update>0 &&  DENSE_THRESHOLD<nrow) {
      nincol=0;
      if (dValue)
	mpt_ft[nincol++]=iPiv;
      for (int i=start;i<end;i++) {
	int iRow=hrowi[i];
	dpermu_ft[iRow]=dluval[i];
	mpt_ft[nincol++]=iRow;
      }
    } else {
      for (int i=start;i<end;i++) {
	int iRow=hrowi[i];
	dpermu_ft[iRow]=dluval[i];
      }
    }
  }
#else
  bool skipStuff = false;
#endif
  if (fact->if_sparse_update>0 &&  DENSE_THRESHOLD<nrow) {
    if (!skipStuff) {
      /* 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 {
    if (!skipStuff) {
      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);
  }
#ifdef CLP_REUSE_ETAS
  fact->nR_etas=save_nR_etas;
  fact->hpivcoR=save_hpivcoR;
  fact->R_etas_start=save_R_etas_start;
#endif
  
  *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 - 2) 
		 + 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);
  }

#ifdef CLP_REUSE_ETAS
  bool skipStuff = (fact->reintro>=0);
  int save_nR_etas=fact->nR_etas;
  int * save_hpivcoR=fact->hpivcoR;
  int * save_R_etas_start=fact->R_etas_start;
  if (skipStuff) {
    // just move
    int * putSeq = fact->xrsadr+2*fact->nrowmx+2;
    int * position = putSeq+fact->maxinv;
    int * putStart = position+fact->maxinv;
    memset(dwork1_ft,0,nincol*sizeof(double));
    int iPiv=fact->reintro;
    int start=putStart[iPiv]&0x7fffffff;
    int end=putStart[iPiv+1]&0x7fffffff;
    double * COIN_RESTRICT dluval	= fact->xeeadr;
    int * COIN_RESTRICT hrowi	= fact->xeradr;
    double dValue;
    if (fact->reintro<fact->nrow) {
      iPiv++;
      dValue=1.0/dluval[start++];
    } else {
      iPiv=hrowi[--end];
      dValue=dluval[end];
      start++;
      int ndoSkip=0;
      for (int i=fact->nrow;i<fact->reintro;i++) {
	if ((putStart[i]&0x80000000)==0)
	ndoSkip++;
      }
      fact->nR_etas-=ndoSkip;
      fact->hpivcoR+=ndoSkip;
      fact->R_etas_start+=ndoSkip;
    }
    dwork1[iPiv]=dValue;
    if (fact->if_sparse_update>0 &&  DENSE_THRESHOLD<nrow) {
      nincol=0;
      if (dValue)
	mpt_ft[nincol++]=iPiv;
      for (int i=start;i<end;i++) {
	int iRow=hrowi[i];
	dwork1[iRow]=dluval[i];
	mpt_ft[nincol++]=iRow;
      }
    } else {
      for (int i=start;i<end;i++) {
	int iRow=hrowi[i];
	dwork1[iRow]=dluval[i];
      }
    }
  }
#else
  bool skipStuff = false;
#endif
  if (fact->if_sparse_update>0 &&  DENSE_THRESHOLD<nrow) {
    if (!skipStuff) {
      /* 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 {
    if (!skipStuff) {
      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);
      }
    }
    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;
    }
  }
#ifdef CLP_REUSE_ETAS
  fact->nR_etas=save_nR_etas;
  fact->hpivcoR=save_hpivcoR;
  fact->R_etas_start=save_R_etas_start;
#endif
  
  
  /* 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);

}