CoinOslC.h

Go to the documentation of this file.

/* $Id: CoinOslC.h 1582 2013-04-06 14:33:07Z stefan $ */
#ifndef COIN_OSL_C_INCLUDE
/*
  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).
*/
#define COIN_OSL_C_INCLUDE

#ifndef CLP_OSL
#define CLP_OSL 0
#endif
#define C_EKK_GO_SPARSE 200

#ifdef HAVE_ENDIAN_H
#include <endian.h>
#if __BYTE_ORDER == __LITTLE_ENDIAN
#define INTEL
#endif
#endif

#include <math.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>

#define SPARSE_UPDATE
#define NO_SHIFT
#include "CoinHelperFunctions.hpp"

#include <stddef.h>
#ifdef __cplusplus
extern "C"{
#endif

int c_ekkbtrn( register const EKKfactinfo *fact,
            double *dwork1,
            int * mpt,int first_nonzero);
int c_ekkbtrn_ipivrw( register const EKKfactinfo *fact,
                   double *dwork1,
                   int * mpt, int ipivrw,int * spare);

int c_ekketsj( register /*const*/ EKKfactinfo *fact,
            double *dwork1,
            int *mpt2, double dalpha, int orig_nincol,
            int npivot, int *nuspikp,
            const int ipivrw, int * spare);
int c_ekkftrn( register const EKKfactinfo *fact, 
            double *dwork1,
            double * dpermu,int * mpt, int numberNonZero);

int c_ekkftrn_ft( register EKKfactinfo *fact, 
               double *dwork1, int *mpt, int *nincolp);
void c_ekkftrn2( register EKKfactinfo *fact, double *dwork1,
              double * dpermu1,int * mpt1, int *nincolp,
             double *dwork1_ft, int *mpt_ft, int *nincolp_ft);

int c_ekklfct( register EKKfactinfo *fact);
int c_ekkslcf( register const EKKfactinfo *fact);
inline void c_ekkscpy(int n, const int *marr1,int *marr2)
{ CoinMemcpyN(marr1,n,marr2);} 
inline void c_ekkdcpy(int n, const double *marr1,double *marr2)
{ CoinMemcpyN(marr1,n,marr2);} 
int c_ekk_IsSet(const int * array,int bit);
void c_ekk_Set(int * array,int bit);
void c_ekk_Unset(int * array,int bit);

void c_ekkzero(int length, int n, void * array);
inline void c_ekkdzero(int n, double *marray)
{CoinZeroN(marray,n);}
inline void c_ekkizero(int n, int *marray)
{CoinZeroN(marray,n);}
inline void c_ekkczero(int n, char *marray)
{CoinZeroN(marray,n);}
#ifdef __cplusplus
          }
#endif
 
#define c_ekkscpy_0_1(s,ival,array) CoinFillN(array,s,ival)
#define c_ekks1cpy( n,marr1,marr2)  CoinMemcpyN(marr1,n, marr2)
void clp_setup_pointers(EKKfactinfo * fact);
void clp_memory(int type);
double * clp_double(int number_entries);
int * clp_int(int number_entries);
void * clp_malloc(int number_entries);
void clp_free(void * oldArray);

#define SLACK_VALUE -1.0
#define C_EKK_REMOVE_LINK(hpiv,hin,link,ipivot) \
  {                                             \
    int ipre = link[ipivot].pre;                \
    int isuc = link[ipivot].suc;                \
    if (ipre > 0) {                             \
      link[ipre].suc = isuc;                    \
    }                                           \
    if (ipre <= 0) {                            \
      hpiv[hin[ipivot]] = isuc;                 \
    }                                           \
    if (isuc > 0) {                             \
      link[isuc].pre = ipre;                    \
    }                                           \
  }

#define C_EKK_ADD_LINK(hpiv,nzi,link, npr)      \
  {                                             \
    int ifiri = hpiv[nzi];                      \
    hpiv[nzi] = npr;                            \
    link[npr].suc = ifiri;                      \
    link[npr].pre = 0;                          \
    if (ifiri != 0) {                           \
      link[ifiri].pre = npr;                    \
    }                                           \
  }
#include <assert.h>
#ifdef  NO_SHIFT

#define SHIFT_INDEX(limit)      (limit)
#define UNSHIFT_INDEX(limit)    (limit)
#define SHIFT_REF(arr,ind)      (arr)[ind]

#else

#define SHIFT_INDEX(limit)      ((limit)<<3)
#define UNSHIFT_INDEX(limit)    ((unsigned int)(limit)>>3)
#define SHIFT_REF(arr,ind)      (*(double*)((char*)(arr) + (ind)))

#endif

#ifdef INTEL
#define NOT_ZERO(x)     (((*((reinterpret_cast<unsigned char *>(&x))+7)) & 0x7F) != 0)
#else
#define NOT_ZERO(x)     ((x) != 0.0)
#endif

#define SWAP(type,_x,_y)        { type _tmp = (_x); (_x) = (_y); (_y) = _tmp;}

#define UNROLL_LOOP_BODY1(code)                 \
  {{code}}
#define UNROLL_LOOP_BODY2(code)                 \
  {{code} {code}}
#define UNROLL_LOOP_BODY4(code)                 \
  {{code} {code} {code} {code}}
#endif
#ifdef COIN_OSL_CMFC
/*     Return codes in IRTCOD/IRTCOD are */
/*     4: numerical problems */
/*     5: not enough space in row file */
/*     6: not enough space in column file */
/*    23: system error at label 320 */
{
#if 1
  int *hcoli    = fact->xecadr;
  double *dluval        = fact->xeeadr;
  double *dvalpv = fact->kw3adr;
  int *mrstrt   = fact->xrsadr;
  int *hrowi    = fact->xeradr;
  int *mcstrt   = fact->xcsadr;
  int *hinrow   = fact->xrnadr;
  int *hincol   = fact->xcnadr;
  int *hpivro   = fact->krpadr; 
  int *hpivco   = fact->kcpadr;
#endif
  int nnentl    = fact->nnentl;
  int nnentu    = fact->nnentu;
  int kmxeta    = fact->kmxeta;
  int xnewro    = *xnewrop;
  int ncompactions      = *ncompactionsp;

  MACTION_T *maction = reinterpret_cast<MACTION_T*>(maction_void);

  int i, j, k;
  double d1;
  int j1, j2;
  int jj, kk, kr, nz, jj1, jj2, kce, kcs, kqq, npr;
  int fill, naft;
  int enpr;
  int nres, npre;
  int knpr, irow, iadd32, ibase;
  double pivot;
  int count, nznpr;
  int nlast, epivr1;
  int kipis;
  double dpivx;
  int kipie, kcpiv, knprs, knpre;
  bool cancel;
  double multip, elemnt;
  int ipivot, jpivot, epivro, epivco, lstart, nfirst;
  int nzpivj, kfill, kstart;
  int nmove, ileft;
#ifndef C_EKKCMFY
  int iput, nspare;
  int noRoomForDense=0;
  int if_sparse_update=fact->if_sparse_update;
  int ifdens = 0;
#endif
  int irtcod    = 0;
  const int nrow        = fact->nrow;

  /* Parameter adjustments */
  --maction;

  /* Function Body */
  lstart = nnetas - nnentl + 1;
  for (i = lstart; i <= nnetas; ++i) {
      hrowi[i] = SHIFT_INDEX(hcoli[i]);
  }

  for (i = 1; i <= nrow; ++i) {
    maction[i] = 0;
    mwork[i].pre = i - 1;
    mwork[i].suc = i + 1;
  }

  iadd32 = 0;
  nlast = nrow;
  nfirst = 1;
  mwork[1].pre = nrow;
  mwork[nrow].suc = 1;

  for (count = 1; count <= nrow; ++count) {

    /* Pick column singletons */
    if (! (hpivco[1] <= 0)) {
      int small_pivot = c_ekkcsin(fact,
                                 rlink, clink,
                                    nsingp);

      if (small_pivot) {
        irtcod = 7; /* pivot too small */
        if (fact->invok >= 0) {
          goto L1050;
        }
      }
      if (fact->npivots >= nrow) {
        goto L1050;
      }
    }

    /* Pick row singletons */
    if (! (hpivro[1] <= 0)) {
      irtcod = c_ekkrsin(fact,
                         rlink, clink,
                         mwork,nfirst,
                         nsingp,
                         
                     &xnewco, &xnewro,
                     &nnentu,
                     &kmxeta, &ncompactions,
                         &nnentl);
        if (irtcod != 0) {
          if (irtcod < 0 || fact->invok >= 0) {
            /* -5 */
            goto L1050;
          }
          /* ASSERT:  irtcod == 7 - pivot too small */
          /* why don't we return with an error? */          
        }
        if (fact->npivots >= nrow) {
            goto L1050;
        }
        lstart = nnetas - nnentl + 1;
    }

    /* Find a pivot element */
    irtcod = c_ekkfpvt(fact,
                      rlink, clink,
                     nsingp, xrejctp, &ipivot, &jpivot);
    if (irtcod != 0) {
      /* irtcod == 10 */
        goto L1050;
    }
    /*        Update list structures and prepare for numerical phase */
    c_ekkprpv(fact, rlink, clink,
                     *xrejctp, ipivot, jpivot);

    epivco = hincol[jpivot];
    ++fact->xnetal;
    mcstrt[fact->xnetal] = lstart - 1;
    hpivco[fact->xnetal] = ipivot;
    epivro = hinrow[ipivot];
    epivr1 = epivro - 1;
    kipis = mrstrt[ipivot];
    pivot = dluval[kipis];
    dpivx = 1. / pivot;
    kipie = kipis + epivr1;
    ++kipis;
#ifndef C_EKKCMFY
    {
      double size = nrow - fact->npivots;
      if (size > GO_DENSE && (nnentu - fact->nuspike) * GO_DENSE_RATIO > size * size) {
        /* say going to dense coding */
        if (*nsingp == 0) {
          ifdens = 1;
        }
      }
    }
#endif
    /* copy the pivot row entries into dvalpv */
    /* the maction array tells us the index into dvalpv for a given row */
    /* the alternative would be using a large array of doubles */
    for (k = kipis; k <= kipie; ++k) {
      irow = hcoli[k];
      dvalpv[k - kipis + 1] = dluval[k];
      maction[irow] = static_cast<MACTION_T>(k - kipis + 1);
    }

    /* Loop over nonzeros in pivot column */
    kcpiv = mcstrt[jpivot] - 1;
    for (nzpivj = 1; nzpivj <= epivco; ++nzpivj) {
      ++kcpiv;
      npr = hrowi[kcpiv];
      hrowi[kcpiv] = 0; /* zero out for possible compaction later on */

      --hincol[jpivot];

      ++mcstrt[jpivot];
      /* loop invariant:  kcpiv == mcstrt[jpivot] - 1 */

      --hinrow[npr];
      enpr = hinrow[npr];
      knprs = mrstrt[npr];
      knpre = knprs + enpr;

      /* Search for element to be eliminated */
      knpr = knprs;
      while (1) {
          UNROLL_LOOP_BODY4({
            if (jpivot == hcoli[knpr]) {
              break;
            }
            knpr++;
          });
      }

      multip = -dluval[knpr] * dpivx;

      /* swap last entry with pivot */
      dluval[knpr] = dluval[knpre];
      hcoli[knpr] = hcoli[knpre];
      --knpre;

#if     1
      /* MONSTER_UNROLLED_CODE - see below */
      kfill = epivr1 - (knpre - knprs + 1);
      nres = ((knpre - knprs + 1) & 1) + knprs;
      cancel = false;
      d1 = 1e33;
      j1 = hcoli[nres];

      if (nres != knprs) {
        j = hcoli[knprs];
        if (maction[j] == 0) {
          ++kfill;
        } else {
          jj = maction[j];
          maction[j] = static_cast<MACTION_T>(-maction[j]);
          dluval[knprs] += multip * dvalpv[jj];
          d1 = fabs(dluval[knprs]);
        }
      }
      j2 = hcoli[nres + 1];
      jj1 = maction[j1];
      for (kr = nres; kr < knpre; kr += 2) {
        jj2 = maction[j2];
        if ( (jj1 == 0)) {
          ++kfill;
        } else {
          maction[j1] = static_cast<MACTION_T>(-maction[j1]);
          dluval[kr] += multip * dvalpv[jj1];
          cancel = cancel || ! (fact->zeroTolerance < d1);
          d1 = fabs(dluval[kr]);
        }
        j1 = hcoli[kr + 2];
        if ( (jj2 == 0)) {
          ++kfill;
        } else {
          maction[j2] = static_cast<MACTION_T>(-maction[j2]);
          dluval[kr + 1] += multip * dvalpv[jj2];
          cancel = cancel || ! (fact->zeroTolerance < d1);
          d1 = fabs(dluval[kr + 1]);
        }
        jj1 = maction[j1];
        j2 = hcoli[kr + 3];
      }
      cancel = cancel || ! (fact->zeroTolerance < d1);
#else
      /*
       * This is apparently what the above code does.
       * In addition to being unrolled, the assignments to j[12] and jj[12]
       * are shifted so that the result of dereferencing maction doesn't
       * have to be used immediately afterwards for the branch test.
       * This would would cause a pipeline delay.  (The apparent dereference
       * of hcoli will be removed by the compiler using strength reduction).
       *
       * loop through the entries in the row being processed,
       * flipping the sign of the maction entries as we go along.
       * Afterwards, we look for positive entries to see what pivot
       * row entries will cause fill-in.  We count the number of fill-ins, too.
       * "cancel" says if the result of combining the pivot row with this one
       * causes an entry to get too small; if so, we discard those entries.
       */
      kfill = epivr1 - (knpre - knprs + 1);
      cancel = false;

      for (kr = knprs; kr <= knpre; kr++) {
        j1 = hcoli[kr];
        jj1 = maction[j1];
        if ( (jj1 == 0)) {
          /* no entry - this pivot column entry will have to be added */
          ++kfill;
        } else {
          /* there is an entry for this column in the pivot row */
          maction[j1] = -maction[j1];
          dluval[kr] += multip * dvalpv[jj1];
          d1 = fabs(dluval[kr]);
          cancel = cancel || ! (fact->zeroTolerance < d1);
        }
      }
#endif
      kstart = knpre;
      fill = kfill;
      
      if (cancel) {
        /* KSTART is used as a stack pointer for nonzeros in factored row */
        kstart = knprs - 1;
        for (kr = knprs; kr <= knpre; ++kr) {
          j = hcoli[kr];
          if (fabs(dluval[kr]) > fact->zeroTolerance) {
            ++kstart;
            dluval[kstart] = dluval[kr];
            hcoli[kstart] = j;
          } else {
            /* Remove element from column file */
            --nnentu;
            --hincol[j];
            --enpr;
            kcs = mcstrt[j];
            kce = kcs + hincol[j];
            for (kk = kcs; kk <= kce; ++kk) {
              if (hrowi[kk] == npr) {
                hrowi[kk] = hrowi[kce];
                hrowi[kce] = 0;
                break;
              }
            }
            /* ASSERT !(kk>kce) */
          }
        }
        knpre = kstart;
      }
      /* Fill contains an upper bound on the amount of fill-in */
      if (fill == 0) {
        for (k = kipis; k <= kipie; ++k) {
          maction[hcoli[k]] = static_cast<MACTION_T>(-maction[hcoli[k]]);
        }
      }
      else {
        naft = mwork[npr].suc;
        kqq = mrstrt[naft] - knpre - 1;
        
        if (fill > kqq) {
          /* Fill-in exceeds space left. Check if there is enough */
          /* space in row file for the new row. */
          nznpr = enpr + fill;
          if (! (xnewro + nznpr + 1 < lstart)) {
            if (! (nnentu + nznpr + 1 < lstart)) {
              irtcod = -5;
              goto L1050;
            }
            /* idea 1 is to compress every time xnewro increases by x thousand */
            /* idea 2 is to copy nucleus rows with a reasonable gap */
            /* then copy each row down when used */
            /* compressions would just be 1 remainder which eventually will */
            /* fit in cache */
            {
              int iput = c_ekkrwcs(fact,dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
              kmxeta += xnewro - iput ;
              xnewro = iput - 1;
              ++ncompactions;
            }
            
            kipis = mrstrt[ipivot] + 1;
            kipie = kipis + epivr1 - 1;
            knprs = mrstrt[npr];
          }
          
          /* I think this assignment should be inside the previous if-stmt */
          /* otherwise, it does nothing */
          /*assert(knpre == knprs + enpr - 1);*/
          knpre = knprs + enpr - 1; 
          
          /*
           * copy this row to the end of the row file and adjust its links.
           * The links keep track of the order of rows in memory.
           * Rows are only moved from the middle all the way to the end.
           */
          if (npr != nlast) {
            npre = mwork[npr].pre;
            if (npr == nfirst) {
              nfirst = naft;
            }
            /*             take out of chain */
            mwork[naft].pre = npre;
            mwork[npre].suc = naft;
            /*             and put in at end */
            mwork[nfirst].pre = npr;
            mwork[nlast].suc = npr;
            mwork[npr].pre = nlast;
            mwork[npr].suc = nfirst;
            nlast = npr;
            kstart = xnewro;
            mrstrt[npr] = kstart + 1;
            nmove = knpre - knprs + 1;
            ibase = kstart + 1 - knprs;
            for (kr = knprs; kr <= knpre; ++kr) {
              dluval[ibase + kr] = dluval[kr];
              hcoli[ibase + kr] = hcoli[kr];
            }
            kstart += nmove;
          } else {
            kstart = knpre;
          }
          
          /* extra space ? */
          /*
           * The mystery of iadd32.
           * This code assigns to xnewro, possibly using iadd32.
           * However, in that case xnewro is assigned to just after
           * the for-loop below, and there is no intervening reference.
           * Therefore, I believe that this code can be entirely eliminated;
           * it is the leftover of an interrupted or dropped experiment.
           * Presumably, this was trying to implement the ideas about
           * padding expressed above.
           */
          if (iadd32 != 0) {
            xnewro += iadd32;
          } else {
            if (kstart + (nrow << 1) + 100 < lstart) {
              ileft = ((nrow - fact->npivots + 32) & -32);
              if (kstart + ileft * ileft + 32 < lstart) {
                iadd32 = ileft;
                xnewro = CoinMax(kstart,xnewro);
                xnewro = (xnewro & -32) + ileft;
              } else {
                xnewro = ((kstart + 31) & -32);
              }
            } else {
              xnewro = kstart;
            }
          }
          
          hinrow[npr] = enpr;
        } else if (! (nnentu + kqq + 2 < lstart)) {
          irtcod = -5;
          goto L1050;
        }
        /* Scan pivot row again to generate fill in. */
        for (kr = kipis; kr <= kipie; ++kr) {
          j = hcoli[kr];
          jj = maction[j];
          if (jj >0) {
            elemnt = multip * dvalpv[jj];
            if (fabs(elemnt) > fact->zeroTolerance) {
              ++kstart;
              dluval[kstart] = elemnt;
              //printf("pivot %d at %d col %d el %g\n",
              // npr,kstart,j,elemnt);
              hcoli[kstart] = j;
              ++nnentu;
              nz = hincol[j];
              kcs = mcstrt[j];
              kce = kcs + nz - 1;
              if (kce == xnewco) {
                if (xnewco + 1 >= lstart) {
                  if (xnewco + nz + 1 >= lstart) {
                    /*                  Compress column file */
                    if (nnentu + nz + 1 < lstart) {
                      xnewco = c_ekkclco(fact,hrowi, mcstrt, hincol, xnewco);
                      ++ncompactions;
                      
                      kcpiv = mcstrt[jpivot] - 1;
                      kcs = mcstrt[j];
                      /*                  HINCOL MAY HAVE CHANGED? (JJHF) ??? */
                      nz = hincol[j];
                      kce = kcs + nz - 1;
                    } else {
                      irtcod = -5;
                      goto L1050;
                    }
                  }
                  /*              Copy column */
                  mcstrt[j] = xnewco + 1;
                  ibase = mcstrt[j] - kcs;
                  for (kk = kcs; kk <= kce; ++kk) {
                    hrowi[ibase + kk] = hrowi[kk];
                    hrowi[kk] = 0;
                  }
                  kce = xnewco + kce - kcs + 1;
                  xnewco = kce + 1;
                } else {
                  ++xnewco;
                }
              } else if (hrowi[kce + 1] != 0) {
                /* here we use the fact that hrowi entries not "in use" are zeroed */
                if (xnewco + nz + 1 >= lstart) {
                  /* Compress column file */
                  if (nnentu + nz + 1 < lstart) {
                    xnewco = c_ekkclco(fact,hrowi, mcstrt, hincol, xnewco);
                    ++ncompactions;
                    
                    kcpiv = mcstrt[jpivot] - 1;
                    kcs = mcstrt[j];
                    /*                  HINCOL MAY HAVE CHANGED? (JJHF) ??? */
                    nz = hincol[j];
                    kce = kcs + nz - 1;
                  } else {
                    irtcod = -5;
                    goto L1050;
                  }
                }
                /* move the column to the end of the column file */
                mcstrt[j] = xnewco + 1;
                ibase = mcstrt[j] - kcs;
                for (kk = kcs; kk <= kce; ++kk) {
                  hrowi[ibase + kk] = hrowi[kk];
                  hrowi[kk] = 0;
                }
                kce = xnewco + kce - kcs + 1;
                xnewco = kce + 1;
              }
              /* store element */
              hrowi[kce + 1] = npr;
              hincol[j] = nz + 1;
            }
          } else {
            maction[j] = static_cast<MACTION_T>(-maction[j]);
          }
        }
        if (fill > kqq) {
          xnewro = kstart;
        }
      }
      hinrow[npr] = kstart - mrstrt[npr] + 1;
      /* Check if row or column file needs compression */
      if (! (xnewco + 1 < lstart)) {
        xnewco = c_ekkclco(fact,hrowi, mcstrt, hincol, xnewco);
        ++ncompactions;
        
        kcpiv = mcstrt[jpivot] - 1;
      }
      if (! (xnewro + 1 < lstart)) {
        int iput = c_ekkrwcs(fact,dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
        kmxeta += xnewro - iput ;
        xnewro = iput - 1;
        ++ncompactions;
        
        kipis = mrstrt[ipivot] + 1;
        kipie = kipis + epivr1 - 1;
      }
      /* Store elementary row transformation */
      ++nnentl;
      --nnentu;
      --lstart;
      dluval[lstart] = multip;
      
      hrowi[lstart] = SHIFT_INDEX(npr);
#define INLINE_AFPV 3
      /* We could do this while computing values but
         it makes it much more complex.  At least we should get
         reasonable cache behavior by doing it each row */
#if INLINE_AFPV
      {
        int j;
        int nel, krs;
        int koff;
        int * index;
        double * els;
        nel = hinrow[npr];
        krs = mrstrt[npr];
        index=&hcoli[krs];
        els=&dluval[krs];
#if INLINE_AFPV<3
#if INLINE_AFPV==1
        double maxaij = 0.0;
        koff = 0;
        j=0;
        while (j<nel) {
          double d = fabs(els[j]);
          if (maxaij < d) {
            maxaij = d;
            koff=j;
          }
          j++;
        }
#else
        assert (nel);
        koff=0;
        double maxaij=fabs(els[0]);
        for (j=1;j<nel;j++) {
          double d = fabs(els[j]);
          if (maxaij < d) {
            maxaij = d;
            koff=j;
          }
        }
#endif
#else
        double maxaij = 0.0;
        koff = 0;
        j=0;
        if ((nel&1)!=0) {
          maxaij=fabs(els[0]);
          j=1;
        }
        
        while (j<nel) {
          UNROLL_LOOP_BODY2({
              double d = fabs(els[j]);
              if (maxaij < d) {
                maxaij = d;
                koff=j;
              }
              j++;
            });
        }
#endif
        SWAP(int, index[koff], index[0]);
        SWAP(double, els[koff], els[0]);
      }
#endif

      {
        int nzi = hinrow[npr];
        if (nzi > 0) {
          C_EKK_ADD_LINK(hpivro, nzi, rlink, npr);
        }
      }
    }

    /* after pivot move biggest to first in each row */
#if INLINE_AFPV==0
    int nn = mcstrt[fact->xnetal] - lstart + 1;
    c_ekkafpv(hrowi+lstart, hcoli, dluval, mrstrt, hinrow, nn);
#endif

    /* Restore work array */
    for (k = kipis; k <= kipie; ++k) {
      maction[hcoli[k]] = 0;
    }

    if (*xrejctp > 0) {
      for (k = kipis; k <= kipie; ++k) {
        int j = hcoli[k];
        int nzj = hincol[j];
        if (! (nzj <= 0) &&
            ! ((clink[j].pre > nrow && nzj != 1))) {
          C_EKK_ADD_LINK(hpivco, nzj, clink, j);
        }
      }
    } else {
      for (k = kipis; k <= kipie; ++k) {
        int j = hcoli[k];
        int nzj = hincol[j];
        if (! (nzj <= 0)) {
          C_EKK_ADD_LINK(hpivco, nzj, clink, j);
        }
      }
    }
    fact->nuspike += hinrow[ipivot];

    /* Go to dense coding if appropriate */
#ifndef C_EKKCMFY
    if (ifdens != 0) {
      int ndense = nrow - fact->npivots;
      if (! (xnewro + ndense * ndense >= lstart)) {

        /* set up sort order in MACTION */
        c_ekkizero( nrow, reinterpret_cast<int *> (maction+1));
        iput = 0;
        for (i = 1; i <= nrow; ++i) {
          if (clink[i].pre >= 0) {
            ++iput;
            maction[i] = static_cast<short int>(iput);
          }
        }
        /* and get number spare needed */
        nspare = 0;
        for (i = 1; i <= nrow; ++i) {
          if (rlink[i].pre >= 0) {
            nspare = nspare + ndense - hinrow[i];
          }
        }
        if (iput != nrow - fact->npivots) {
          /* must be singular */
          c_ekkizero( nrow, reinterpret_cast<int *> (maction+1));
        } else {
          /* pack down then back up */
          int iput = c_ekkrwcs(fact,dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
          kmxeta += xnewro - iput ;
          xnewro = iput - 1;
          ++ncompactions;
          
          --ncompactions;
          if (xnewro + nspare + ndense * ndense >= lstart) {
            c_ekkizero( nrow, reinterpret_cast<int *> (maction+1));
          }
          else {
            xnewro += nspare;
            c_ekkrwct(fact,dluval, hcoli, mrstrt, hinrow, mwork,
                    rlink, maction, dvalpv,
                    nlast,  xnewro);
            kmxeta += xnewro ;
            if (nnentu + nnentl > nrow * 5 &&
                (ndense*ndense)>(nnentu+nnentl)>>2 &&
                !if_sparse_update) {
              fact->ndenuc = ndense;
            }
            irtcod = c_ekkcmfd(fact,
                             (reinterpret_cast<int*>(dvalpv)+1),
                             rlink, clink,
                             (reinterpret_cast<int*>(maction+1))+1,
                             nnetas,
                             &nnentl, &nnentu,
                             nsingp);
            /* irtcod == 0 || irtcod == 10 */
            /* 10 == found 0.0 pivot */
            goto L1050;
          }
        }
      } else {
        /* say not enough room */
        /*printf("no room %d\n",ndense);*/
        if (1) {
          /* return and increase size of etas if possible */
          if (!noRoomForDense) {
            int etasize =CoinMax(4*fact->nnentu+(nnetas-fact->nnentl)+1000,fact->eta_size);
            noRoomForDense=ndense;
            fact->eta_size=CoinMin(static_cast<int>(1.2*fact->eta_size),etasize);
            if (fact->maxNNetas>0&&fact->eta_size>
                fact->maxNNetas) {
              fact->eta_size=fact->maxNNetas;
            }
          }
        }
      }
    }
#endif  /* C_EKKCMFY */
  }

 L1050:
  {
    int iput = c_ekkrwcs(fact,dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
    kmxeta += xnewro - iput;
    xnewro = iput - 1;
    ++ncompactions;
  }

  nnentu = xnewro;
  /* save order of row copy for c_ekkshfv */
  mwork[nrow+1].pre = nfirst;
  mwork[nrow+1].suc = nlast;

  fact->nnentl = nnentl;
  fact->nnentu = nnentu;
  fact->kmxeta = kmxeta;
  *xnewrop = xnewro;
  *ncompactionsp = ncompactions;

  return (irtcod);
} /* c_ekkcmfc */
#endif