CoinUtils/src/CoinSimpFactorization.cpp

/* $Id: CoinSimpFactorization.cpp 2083 2019-01-06 19:38:09Z unxusr $ */
// Copyright (C) 2008, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).
#if COIN_BIG_INDEX == 0
#include "CoinUtilsConfig.h"

#include <cassert>
#include "CoinPragma.hpp"
#include "CoinSimpFactorization.hpp"
#include "CoinIndexedVector.hpp"
#include "CoinHelperFunctions.hpp"
#include "CoinPackedMatrix.hpp"
#include "CoinFinite.hpp"
#include <stdio.h>

#define ARRAY 0

FactorPointers::FactorPointers(int numRows, int numColumns,
  int *UrowLengths_,
  int *UcolLengths_)
{

  rowMax = new double[numRows];
  double *current = rowMax;
  const double *end = current + numRows;
  for (; current != end; ++current)
    *current = -1.0;

  firstRowKnonzeros = new int[numRows + 1];
  CoinFillN(firstRowKnonzeros, numRows + 1, -1);

  prevRow = new int[numRows];
  nextRow = new int[numRows];
  firstColKnonzeros = new int[numRows + 1];
  memset(firstColKnonzeros, -1, (numRows + 1) * sizeof(int));

  prevColumn = new int[numColumns];
  nextColumn = new int[numColumns];
  newCols = new int[numRows];

  for (int i = numRows - 1; i >= 0; --i) {
    int length = UrowLengths_[i];
    prevRow[i] = -1;
    nextRow[i] = firstRowKnonzeros[length];
    if (nextRow[i] != -1)
      prevRow[nextRow[i]] = i;
    firstRowKnonzeros[length] = i;
  }
  for (int i = numColumns - 1; i >= 0; --i) {
    int length = UcolLengths_[i];
    prevColumn[i] = -1;
    nextColumn[i] = firstColKnonzeros[length];
    if (nextColumn[i] != -1)
      prevColumn[nextColumn[i]] = i;
    firstColKnonzeros[length] = i;
  }
}
FactorPointers::~FactorPointers()
{
  delete[] rowMax;
  delete[] firstRowKnonzeros;
  delete[] prevRow;
  delete[] nextRow;
  delete[] firstColKnonzeros;
  delete[] prevColumn;
  delete[] nextColumn;
  delete[] newCols;
}

//:class CoinSimpFactorization.  Deals with Factorization and Updates
//  CoinSimpFactorization.  Constructor
CoinSimpFactorization::CoinSimpFactorization()
  : CoinOtherFactorization()
{
  gutsOfInitialize();
}

/// Copy constructor
CoinSimpFactorization::CoinSimpFactorization(const CoinSimpFactorization &other)
  : CoinOtherFactorization(other)
{
  gutsOfInitialize();
  gutsOfCopy(other);
}
// Clone
CoinOtherFactorization *
CoinSimpFactorization::clone() const
{
  return new CoinSimpFactorization(*this);
}
/// The real work of constructors etc
void CoinSimpFactorization::gutsOfDestructor()
{
  delete[] elements_;
  delete[] pivotRow_;
  delete[] workArea_;
  elements_ = NULL;
  pivotRow_ = NULL;
  workArea_ = NULL;
  numberRows_ = 0;
  numberColumns_ = 0;
  numberGoodU_ = 0;
  status_ = -1;
  maximumRows_ = 0;
  maximumSpace_ = 0;
  numberSlacks_ = 0;
  firstNumberSlacks_ = 0;
  //
  delete[] denseVector_;
  delete[] workArea2_;
  delete[] workArea3_;
  delete[] vecLabels_;
  delete[] indVector_;

  delete[] auxVector_;
  delete[] auxInd_;

  delete[] vecKeep_;
  delete[] indKeep_;

  delete[] LrowStarts_;
  delete[] LrowLengths_;
  delete[] Lrows_;
  delete[] LrowInd_;

  delete[] LcolStarts_;
  delete[] LcolLengths_;
  delete[] Lcolumns_;
  delete[] LcolInd_;

  delete[] UrowStarts_;
  delete[] UrowLengths_;
#ifdef COIN_SIMP_CAPACITY
  delete[] UrowCapacities_;
#endif
  delete[] Urows_;
  delete[] UrowInd_;

  delete[] prevRowInU_;
  delete[] nextRowInU_;

  delete[] UcolStarts_;
  delete[] UcolLengths_;
#ifdef COIN_SIMP_CAPACITY
  delete[] UcolCapacities_;
#endif
  delete[] Ucolumns_;
  delete[] UcolInd_;
  delete[] prevColInU_;
  delete[] nextColInU_;
  delete[] colSlack_;

  delete[] invOfPivots_;

  delete[] colOfU_;
  delete[] colPosition_;
  delete[] rowOfU_;
  delete[] rowPosition_;
  delete[] secRowOfU_;
  delete[] secRowPosition_;

  delete[] EtaPosition_;
  delete[] EtaStarts_;
  delete[] EtaLengths_;
  delete[] EtaInd_;
  delete[] Eta_;

  denseVector_ = NULL;
  workArea2_ = NULL;
  workArea3_ = NULL;
  vecLabels_ = NULL;
  indVector_ = NULL;

  auxVector_ = NULL;
  auxInd_ = NULL;

  vecKeep_ = NULL;
  indKeep_ = NULL;

  LrowStarts_ = NULL;
  LrowLengths_ = NULL;
  Lrows_ = NULL;
  LrowInd_ = NULL;

  LcolStarts_ = NULL;
  LcolLengths_ = NULL;
  Lcolumns_ = NULL;
  LcolInd_ = NULL;

  UrowStarts_ = NULL;
  UrowLengths_ = NULL;
#ifdef COIN_SIMP_CAPACITY
  UrowCapacities_ = NULL;
#endif
  Urows_ = NULL;
  UrowInd_ = NULL;

  prevRowInU_ = NULL;
  nextRowInU_ = NULL;

  UcolStarts_ = NULL;
  UcolLengths_ = NULL;
#ifdef COIN_SIMP_CAPACITY
  UcolCapacities_ = NULL;
#endif
  Ucolumns_ = NULL;
  UcolInd_ = NULL;
  prevColInU_ = NULL;
  nextColInU_ = NULL;
  colSlack_ = NULL;

  invOfPivots_ = NULL;

  colOfU_ = NULL;
  colPosition_ = NULL;
  rowOfU_ = NULL;
  rowPosition_ = NULL;
  secRowOfU_ = NULL;
  secRowPosition_ = NULL;

  EtaPosition_ = NULL;
  EtaStarts_ = NULL;
  EtaLengths_ = NULL;
  EtaInd_ = NULL;
  Eta_ = NULL;
}
void CoinSimpFactorization::gutsOfInitialize()
{
  pivotTolerance_ = 1.0e-1;
  zeroTolerance_ = 1.0e-13;
#ifndef COIN_FAST_CODE
  slackValue_ = -1.0;
#endif
  maximumPivots_ = 200;
  relaxCheck_ = 1.0;
  numberRows_ = 0;
  numberColumns_ = 0;
  numberGoodU_ = 0;
  status_ = -1;
  numberPivots_ = 0;
  maximumRows_ = 0;
  maximumSpace_ = 0;
  numberSlacks_ = 0;
  firstNumberSlacks_ = 0;
  elements_ = NULL;
  pivotRow_ = NULL;
  workArea_ = NULL;

  denseVector_ = NULL;
  workArea2_ = NULL;
  workArea3_ = NULL;
  vecLabels_ = NULL;
  indVector_ = NULL;

  auxVector_ = NULL;
  auxInd_ = NULL;

  vecKeep_ = NULL;
  indKeep_ = NULL;

  LrowStarts_ = NULL;
  LrowLengths_ = NULL;
  Lrows_ = NULL;
  LrowInd_ = NULL;

  LcolStarts_ = NULL;
  LcolLengths_ = NULL;
  Lcolumns_ = NULL;
  LcolInd_ = NULL;

  UrowStarts_ = NULL;
  UrowLengths_ = NULL;
#ifdef COIN_SIMP_CAPACITY
  UrowCapacities_ = NULL;
#endif
  Urows_ = NULL;
  UrowInd_ = NULL;

  prevRowInU_ = NULL;
  nextRowInU_ = NULL;

  UcolStarts_ = NULL;
  UcolLengths_ = NULL;
#ifdef COIN_SIMP_CAPACITY
  UcolCapacities_ = NULL;
#endif
  Ucolumns_ = NULL;
  UcolInd_ = NULL;
  prevColInU_ = NULL;
  nextColInU_ = NULL;
  colSlack_ = NULL;

  invOfPivots_ = NULL;

  colOfU_ = NULL;
  colPosition_ = NULL;
  rowOfU_ = NULL;
  rowPosition_ = NULL;
  secRowOfU_ = NULL;
  secRowPosition_ = NULL;

  EtaPosition_ = NULL;
  EtaStarts_ = NULL;
  EtaLengths_ = NULL;
  EtaInd_ = NULL;
  Eta_ = NULL;
}
void CoinSimpFactorization::initialSomeNumbers()
{
  keepSize_ = -1;
  LrowSize_ = -1;
  // LrowCap_ in allocateSomeArrays
  LcolSize_ = -1;
  // LcolCap_ in allocateSomeArrays
  // UrowMaxCap_ in allocateSomeArrays
  UrowEnd_ = -1;
  firstRowInU_ = -1;
  lastRowInU_ = -1;
  firstColInU_ = -1;
  lastColInU_ = -1;
  //UcolMaxCap_ in allocateSomeArrays
  UcolEnd_ = -1;

  EtaSize_ = 0;
  lastEtaRow_ = -1;
  //maxEtaRows_ in allocateSomeArrays
  //EtaMaxCap_ in allocateSomeArrays

  // minIncrease_ in allocateSomeArrays
  updateTol_ = 1.0e12;

  doSuhlHeuristic_ = true;
  maxU_ = -1.0;
  maxGrowth_ = 1.e12;
  maxA_ = -1.0;
  pivotCandLimit_ = 4;
  minIncrease_ = 10;
}

//  ~CoinSimpFactorization.  Destructor
CoinSimpFactorization::~CoinSimpFactorization()
{
  gutsOfDestructor();
}
//  =
CoinSimpFactorization &CoinSimpFactorization::operator=(const CoinSimpFactorization &other)
{
  if (this != &other) {
    gutsOfDestructor();
    gutsOfInitialize();
    gutsOfCopy(other);
  }
  return *this;
}
void CoinSimpFactorization::gutsOfCopy(const CoinSimpFactorization &other)
{
  pivotTolerance_ = other.pivotTolerance_;
  zeroTolerance_ = other.zeroTolerance_;
#ifndef COIN_FAST_CODE
  slackValue_ = other.slackValue_;
#endif
  relaxCheck_ = other.relaxCheck_;
  numberRows_ = other.numberRows_;
  numberColumns_ = other.numberColumns_;
  maximumRows_ = other.maximumRows_;
  maximumSpace_ = other.maximumSpace_;
  numberGoodU_ = other.numberGoodU_;
  maximumPivots_ = other.maximumPivots_;
  numberPivots_ = other.numberPivots_;
  factorElements_ = other.factorElements_;
  status_ = other.status_;
  numberSlacks_ = other.numberSlacks_;
  firstNumberSlacks_ = other.firstNumberSlacks_;
  if (other.pivotRow_) {
    pivotRow_ = new int[2 * maximumRows_ + maximumPivots_];
    memcpy(pivotRow_, other.pivotRow_, (2 * maximumRows_ + numberPivots_) * sizeof(int));
    elements_ = new CoinFactorizationDouble[maximumSpace_];
    memcpy(elements_, other.elements_, (maximumRows_ + numberPivots_) * maximumRows_ * sizeof(CoinFactorizationDouble));
    workArea_ = new CoinFactorizationDouble[maximumRows_];
  } else {
    elements_ = NULL;
    pivotRow_ = NULL;
    workArea_ = NULL;
  }

  keepSize_ = other.keepSize_;

  LrowSize_ = other.LrowSize_;
  LrowCap_ = other.LrowCap_;

  LcolSize_ = other.LcolSize_;
  LcolCap_ = other.LcolCap_;

  UrowMaxCap_ = other.UrowMaxCap_;
  UrowEnd_ = other.UrowEnd_;
  firstRowInU_ = other.firstRowInU_;
  lastRowInU_ = other.lastRowInU_;

  firstColInU_ = other.firstColInU_;
  lastColInU_ = other.lastColInU_;
  UcolMaxCap_ = other.UcolMaxCap_;
  UcolEnd_ = other.UcolEnd_;

  EtaSize_ = other.EtaSize_;
  lastEtaRow_ = other.lastEtaRow_;
  maxEtaRows_ = other.maxEtaRows_;
  EtaMaxCap_ = other.EtaMaxCap_;

  minIncrease_ = other.minIncrease_;
  updateTol_ = other.updateTol_;

  if (other.denseVector_) {
    denseVector_ = new double[maximumRows_];
    memcpy(denseVector_, other.denseVector_, maximumRows_ * sizeof(double));
  } else
    denseVector_ = NULL;
  if (other.workArea2_) {
    workArea2_ = new double[maximumRows_];
    memcpy(workArea2_, other.workArea2_, maximumRows_ * sizeof(double));
  } else
    workArea2_ = NULL;
  if (other.workArea3_) {
    workArea3_ = new double[maximumRows_];
    memcpy(workArea3_, other.workArea3_, maximumRows_ * sizeof(double));
  } else
    workArea3_ = NULL;
  if (other.vecLabels_) {
    vecLabels_ = new int[maximumRows_];
    memcpy(vecLabels_, other.vecLabels_, maximumRows_ * sizeof(int));
  } else
    vecLabels_ = NULL;
  if (other.indVector_) {
    indVector_ = new int[maximumRows_];
    memcpy(indVector_, other.indVector_, maximumRows_ * sizeof(int));
  } else
    indVector_ = NULL;

  if (other.auxVector_) {
    auxVector_ = new double[maximumRows_];
    memcpy(auxVector_, other.auxVector_, maximumRows_ * sizeof(double));
  } else
    auxVector_ = NULL;
  if (other.auxInd_) {
    auxInd_ = new int[maximumRows_];
    memcpy(auxInd_, other.auxInd_, maximumRows_ * sizeof(int));
  } else
    auxInd_ = NULL;

  if (other.vecKeep_) {
    vecKeep_ = new double[maximumRows_];
    memcpy(vecKeep_, other.vecKeep_, maximumRows_ * sizeof(double));
  } else
    vecKeep_ = NULL;
  if (other.indKeep_) {
    indKeep_ = new int[maximumRows_];
    memcpy(indKeep_, other.indKeep_, maximumRows_ * sizeof(int));
  } else
    indKeep_ = NULL;
  if (other.LrowStarts_) {
    LrowStarts_ = new int[maximumRows_];
    memcpy(LrowStarts_, other.LrowStarts_, maximumRows_ * sizeof(int));
  } else
    LrowStarts_ = NULL;
  if (other.LrowLengths_) {
    LrowLengths_ = new int[maximumRows_];
    memcpy(LrowLengths_, other.LrowLengths_, maximumRows_ * sizeof(int));
  } else
    LrowLengths_ = NULL;
  if (other.Lrows_) {
    Lrows_ = new double[other.LrowCap_];
    memcpy(Lrows_, other.Lrows_, other.LrowCap_ * sizeof(double));
  } else
    Lrows_ = NULL;
  if (other.LrowInd_) {
    LrowInd_ = new int[other.LrowCap_];
    memcpy(LrowInd_, other.LrowInd_, other.LrowCap_ * sizeof(int));
  } else
    LrowInd_ = NULL;

  if (other.LcolStarts_) {
    LcolStarts_ = new int[maximumRows_];
    memcpy(LcolStarts_, other.LcolStarts_, maximumRows_ * sizeof(int));
  } else
    LcolStarts_ = NULL;
  if (other.LcolLengths_) {
    LcolLengths_ = new int[maximumRows_];
    memcpy(LcolLengths_, other.LcolLengths_, maximumRows_ * sizeof(int));
  } else
    LcolLengths_ = NULL;
  if (other.Lcolumns_) {
    Lcolumns_ = new double[other.LcolCap_];
    memcpy(Lcolumns_, other.Lcolumns_, other.LcolCap_ * sizeof(double));
  } else
    Lcolumns_ = NULL;
  if (other.LcolInd_) {
    LcolInd_ = new int[other.LcolCap_];
    memcpy(LcolInd_, other.LcolInd_, other.LcolCap_ * sizeof(int));
  } else
    LcolInd_ = NULL;

  if (other.UrowStarts_) {
    UrowStarts_ = new int[maximumRows_];
    memcpy(UrowStarts_, other.UrowStarts_, maximumRows_ * sizeof(int));
  } else
    UrowStarts_ = NULL;
  if (other.UrowLengths_) {
    UrowLengths_ = new int[maximumRows_];
    memcpy(UrowLengths_, other.UrowLengths_, maximumRows_ * sizeof(int));
  } else
    UrowLengths_ = NULL;
#ifdef COIN_SIMP_CAPACITY
  if (other.UrowCapacities_) {
    UrowCapacities_ = new int[maximumRows_];
    memcpy(UrowCapacities_, other.UrowCapacities_, maximumRows_ * sizeof(int));
  } else
    UrowCapacities_ = NULL;
#endif
  if (other.Urows_) {
    Urows_ = new double[other.UrowMaxCap_];
    memcpy(Urows_, other.Urows_, other.UrowMaxCap_ * sizeof(double));
  } else
    Urows_ = NULL;
  if (other.UrowInd_) {
    UrowInd_ = new int[other.UrowMaxCap_];
    memcpy(UrowInd_, other.UrowInd_, other.UrowMaxCap_ * sizeof(int));
  } else
    UrowInd_ = NULL;

  if (other.prevRowInU_) {
    prevRowInU_ = new int[maximumRows_];
    memcpy(prevRowInU_, other.prevRowInU_, maximumRows_ * sizeof(int));
  } else
    prevRowInU_ = NULL;
  if (other.nextRowInU_) {
    nextRowInU_ = new int[maximumRows_];
    memcpy(nextRowInU_, other.nextRowInU_, maximumRows_ * sizeof(int));
  } else
    nextRowInU_ = NULL;

  if (other.UcolStarts_) {
    UcolStarts_ = new int[maximumRows_];
    memcpy(UcolStarts_, other.UcolStarts_, maximumRows_ * sizeof(int));
  } else
    UcolStarts_ = NULL;
  if (other.UcolLengths_) {
    UcolLengths_ = new int[maximumRows_];
    memcpy(UcolLengths_, other.UcolLengths_, maximumRows_ * sizeof(int));
  } else
    UcolLengths_ = NULL;
#ifdef COIN_SIMP_CAPACITY
  if (other.UcolCapacities_) {
    UcolCapacities_ = new int[maximumRows_];
    memcpy(UcolCapacities_, other.UcolCapacities_, maximumRows_ * sizeof(int));
  } else
    UcolCapacities_ = NULL;
#endif
  if (other.Ucolumns_) {
    Ucolumns_ = new double[other.UcolMaxCap_];
    memcpy(Ucolumns_, other.Ucolumns_, other.UcolMaxCap_ * sizeof(double));
  } else
    Ucolumns_ = NULL;
  if (other.UcolInd_) {
    UcolInd_ = new int[other.UcolMaxCap_];
    memcpy(UcolInd_, other.UcolInd_, other.UcolMaxCap_ * sizeof(int));
  } else
    UcolInd_ = NULL;
  if (other.prevColInU_) {
    prevColInU_ = new int[maximumRows_];
    memcpy(prevColInU_, other.prevColInU_, maximumRows_ * sizeof(int));
  } else
    prevColInU_ = NULL;
  if (other.nextColInU_) {
    nextColInU_ = new int[maximumRows_];
    memcpy(nextColInU_, other.nextColInU_, maximumRows_ * sizeof(int));
  } else
    nextColInU_ = NULL;
  if (other.colSlack_) {
    colSlack_ = new int[maximumRows_];
    memcpy(colSlack_, other.colSlack_, maximumRows_ * sizeof(int));
  }

  if (other.invOfPivots_) {
    invOfPivots_ = new double[maximumRows_];
    memcpy(invOfPivots_, other.invOfPivots_, maximumRows_ * sizeof(double));
  } else
    invOfPivots_ = NULL;

  if (other.colOfU_) {
    colOfU_ = new int[maximumRows_];
    memcpy(colOfU_, other.colOfU_, maximumRows_ * sizeof(int));
  } else
    colOfU_ = NULL;
  if (other.colPosition_) {
    colPosition_ = new int[maximumRows_];
    memcpy(colPosition_, other.colPosition_, maximumRows_ * sizeof(int));
  } else
    colPosition_ = NULL;
  if (other.rowOfU_) {
    rowOfU_ = new int[maximumRows_];
    memcpy(rowOfU_, other.rowOfU_, maximumRows_ * sizeof(int));
  } else
    rowOfU_ = NULL;
  if (other.rowPosition_) {
    rowPosition_ = new int[maximumRows_];
    memcpy(rowPosition_, other.rowPosition_, maximumRows_ * sizeof(int));
  } else
    rowPosition_ = NULL;
  if (other.secRowOfU_) {
    secRowOfU_ = new int[maximumRows_];
    memcpy(secRowOfU_, other.secRowOfU_, maximumRows_ * sizeof(int));
  } else
    secRowOfU_ = NULL;
  if (other.secRowPosition_) {
    secRowPosition_ = new int[maximumRows_];
    memcpy(secRowPosition_, other.secRowPosition_, maximumRows_ * sizeof(int));
  } else
    secRowPosition_ = NULL;

  if (other.EtaPosition_) {
    EtaPosition_ = new int[other.maxEtaRows_];
    memcpy(EtaPosition_, other.EtaPosition_, other.maxEtaRows_ * sizeof(int));
  } else
    EtaPosition_ = NULL;
  if (other.EtaStarts_) {
    EtaStarts_ = new int[other.maxEtaRows_];
    memcpy(EtaStarts_, other.EtaStarts_, other.maxEtaRows_ * sizeof(int));
  } else
    EtaStarts_ = NULL;
  if (other.EtaLengths_) {
    EtaLengths_ = new int[other.maxEtaRows_];
    memcpy(EtaLengths_, other.EtaLengths_, other.maxEtaRows_ * sizeof(int));
  } else
    EtaLengths_ = NULL;
  if (other.EtaInd_) {
    EtaInd_ = new int[other.EtaMaxCap_];
    memcpy(EtaInd_, other.EtaInd_, other.EtaMaxCap_ * sizeof(int));
  } else
    EtaInd_ = NULL;
  if (other.Eta_) {
    Eta_ = new double[other.EtaMaxCap_];
    memcpy(Eta_, other.Eta_, other.EtaMaxCap_ * sizeof(double));
  } else
    Eta_ = NULL;

  doSuhlHeuristic_ = other.doSuhlHeuristic_;
  maxU_ = other.maxU_;
  maxGrowth_ = other.maxGrowth_;
  maxA_ = other.maxA_;
  pivotCandLimit_ = other.pivotCandLimit_;
}

//  getAreas.  Gets space for a factorization
//called by constructors
void CoinSimpFactorization::getAreas(int numberOfRows,
  int numberOfColumns,
  int,
  int)
{

  numberRows_ = numberOfRows;
  numberColumns_ = numberOfColumns;
  int size = numberRows_ * (numberRows_ + CoinMax(maximumPivots_, (numberRows_ + 1) >> 1));
  if (size > maximumSpace_) {
    delete[] elements_;
    elements_ = new CoinFactorizationDouble[size];
    maximumSpace_ = size;
  }
  if (numberRows_ > maximumRows_) {
    maximumRows_ = numberRows_;
    delete[] pivotRow_;
    delete[] workArea_;
    pivotRow_ = new int[2 * maximumRows_ + maximumPivots_];
    workArea_ = new CoinFactorizationDouble[maximumRows_];
    allocateSomeArrays();
  }
}

//  preProcess.
void CoinSimpFactorization::preProcess()
{
  int put = numberRows_ * numberRows_;
  int *indexRow = reinterpret_cast< int * >(elements_ + put);
  int *starts = reinterpret_cast< int * >(pivotRow_);
  initialSomeNumbers();

  // compute sizes for Urows_ and Ucolumns_
  //for ( int row=0; row < numberRows_; ++row )
  //UrowLengths_[row]=0;
  int k = 0;
  for (int column = 0; column < numberColumns_; ++column) {
    UcolStarts_[column] = k;
    //for (int j=starts[column];j<starts[column+1];j++) {
    //  int iRow = indexRow[j];
    //  ++UrowLengths_[iRow];
    //	}
    UcolLengths_[column] = starts[column + 1] - starts[column];
#ifdef COIN_SIMP_CAPACITY
    UcolCapacities_[column] = numberRows_;
#endif
    //k+=UcolLengths_[column]+minIncrease_;
    k += numberRows_;
  }

  // create space for Urows_
  k = 0;
  for (int row = 0; row < numberRows_; ++row) {
    prevRowInU_[row] = row - 1;
    nextRowInU_[row] = row + 1;
    UrowStarts_[row] = k;
    //k+=UrowLengths_[row]+minIncrease_;
    k += numberRows_;
    UrowLengths_[row] = 0;
#ifdef COIN_SIMP_CAPACITY
    UrowCapacities_[row] = numberRows_;
#endif
  }
  UrowEnd_ = k;
  nextRowInU_[numberRows_ - 1] = -1;
  firstRowInU_ = 0;
  lastRowInU_ = numberRows_ - 1;
  maxA_ = -1.0;
  // build Ucolumns_ and Urows_
  int colBeg, colEnd;
  for (int column = 0; column < numberColumns_; ++column) {
    prevColInU_[column] = column - 1;
    nextColInU_[column] = column + 1;
    k = 0;
    colBeg = starts[column];
    colEnd = starts[column + 1];
    // identify slacks
    if (colEnd == colBeg + 1 && elements_[colBeg] == slackValue_)
      colSlack_[column] = 1;
    else
      colSlack_[column] = 0;
    //
    for (int j = colBeg; j < colEnd; j++) {
      // Ucolumns_
      int iRow = indexRow[j];
      UcolInd_[UcolStarts_[column] + k] = iRow;
      ++k;
      // Urow_
      int ind = UrowStarts_[iRow] + UrowLengths_[iRow];
      UrowInd_[ind] = column;
      Urows_[ind] = elements_[j];
      //maxA_=CoinMax( maxA_, fabs(Urows_[ind]) );
      ++UrowLengths_[iRow];
    }
  }
  nextColInU_[numberColumns_ - 1] = -1;
  firstColInU_ = 0;
  lastColInU_ = numberColumns_ - 1;

  // Initialize L
  LcolSize_ = 0;
  //LcolCap_=numberRows_*minIncrease_;
  memset(LrowStarts_, -1, numberRows_ * sizeof(int));
  memset(LrowLengths_, 0, numberRows_ * sizeof(int));
  memset(LcolStarts_, -1, numberRows_ * sizeof(int));
  memset(LcolLengths_, 0, numberRows_ * sizeof(int));

  // initialize permutations
  for (int i = 0; i < numberRows_; ++i) {
    rowOfU_[i] = i;
    rowPosition_[i] = i;
  }
  for (int i = 0; i < numberColumns_; ++i) {
    colOfU_[i] = i;
    colPosition_[i] = i;
  }
  //

  doSuhlHeuristic_ = true;
}

void CoinSimpFactorization::factorize(int numberOfRows,
  int numberOfColumns,
  const int colStarts[],
  const int indicesRow[],
  const double elements[])
{
  int maximumL = 0;
  int maximumU = 0;
  getAreas(numberOfRows, numberOfColumns, maximumL, maximumU);

  int put = numberRows_ * numberRows_;
  int *indexRow = reinterpret_cast< int * >(elements_ + put);
  int *starts = reinterpret_cast< int * >(pivotRow_);
  for (int column = 0; column <= numberColumns_; ++column) {
    starts[column] = colStarts[column];
  }
  const int limit = colStarts[numberColumns_];
  for (int i = 0; i < limit; ++i) {
    indexRow[i] = indicesRow[i];
    elements_[i] = elements[i];
  }

  preProcess();
  factor();
}

//Does factorization
int CoinSimpFactorization::factor()
{
  numberPivots_ = 0;
  status_ = 0;

  FactorPointers pointers(numberRows_, numberColumns_, UrowLengths_, UcolLengths_);
  int rc = mainLoopFactor(pointers);
  // rc=0 success
  if (rc != 0) { // failure
    status_ = -1;
    //return status_; // failure
  }
  //if ( rc == -3 ) {  // numerical instability
  //    status_ = -1;
  //    return status_;
  //}

  //copyLbyRows();
  copyUbyColumns();
  copyRowPermutations();
  firstNumberSlacks_ = numberSlacks_;
  // row permutations
  if (status_ == -1 || numberColumns_ < numberRows_) {
    for (int j = 0; j < numberRows_; j++)
      pivotRow_[j + numberRows_] = rowOfU_[j];
    for (int j = 0; j < numberRows_; j++) {
      int k = pivotRow_[j + numberRows_];
      pivotRow_[k] = j;
    }
  } else // no permutations
    for (int j = 0; j < numberRows_; j++) {
      pivotRow_[j] = j;
      pivotRow_[j + numberRows_] = j;
    }

  return status_;
}
//

// Makes a non-singular basis by replacing variables
void CoinSimpFactorization::makeNonSingular(int *sequence, int numberColumns)
{
  // Replace bad ones by correct slack
  int *workArea = reinterpret_cast< int * >(workArea_);
  int i;
  for (i = 0; i < numberRows_; i++)
    workArea[i] = -1;
  for (i = 0; i < numberGoodU_; i++) {
    int iOriginal = pivotRow_[i + numberRows_];
    workArea[iOriginal] = i;
    //workArea[i]=iOriginal;
  }
  int lastRow = -1;
  for (i = 0; i < numberRows_; i++) {
    if (workArea[i] == -1) {
      lastRow = i;
      break;
    }
  }
  assert(lastRow >= 0);
  for (i = numberGoodU_; i < numberRows_; i++) {
    assert(lastRow < numberRows_);
    // Put slack in basis
    sequence[i] = lastRow + numberColumns;
    lastRow++;
    for (; lastRow < numberRows_; lastRow++) {
      if (workArea[lastRow] == -1)
        break;
    }
  }
}
// Does post processing on valid factorization - putting variables on correct rows
void CoinSimpFactorization::postProcess(const int *sequence, int *pivotVariable)
{
  for (int i = 0; i < numberRows_; i++) {
    int k = sequence[i];
    pivotVariable[pivotRow_[i + numberRows_]] = k;
  }
}
/* Replaces one Column to basis,
   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
   If checkBeforeModifying is true will do all accuracy checks
   before modifying factorization.  Whether to set this depends on
   speed considerations.  You could just do this on first iteration
   after factorization and thereafter re-factorize
   partial update already in U */
int CoinSimpFactorization::replaceColumn(CoinIndexedVector *,
  int pivotRow,
  double pivotCheck,
  bool,
  double)
{
  if (numberPivots_ == maximumPivots_)
    return 3;

  double pivotValue = pivotCheck;
  if (fabs(pivotValue) < zeroTolerance_)
    return 2;
  int realPivotRow = pivotRow_[pivotRow];

  LUupdate(pivotRow);

  pivotRow_[2 * numberRows_ + numberPivots_] = realPivotRow;
  numberPivots_++;

  return 0;
}
int CoinSimpFactorization::updateColumn(CoinIndexedVector *regionSparse,
  CoinIndexedVector *regionSparse2,
  bool noPermute) const
{
  return upColumn(regionSparse, regionSparse2, noPermute, false);
}
int CoinSimpFactorization::updateColumnFT(CoinIndexedVector *regionSparse,
  CoinIndexedVector *regionSparse2,
  bool noPermute)
{

  int rc = upColumn(regionSparse, regionSparse2, noPermute, true);
  return rc;
}

int CoinSimpFactorization::upColumn(CoinIndexedVector *regionSparse,
  CoinIndexedVector *regionSparse2,
  bool, bool save) const
{
  assert(numberRows_ == numberColumns_);
  double *region2 = regionSparse2->denseVector();
  int *regionIndex = regionSparse2->getIndices();
  int numberNonZero = regionSparse2->getNumElements();
  double *region = regionSparse->denseVector();
  if (!regionSparse2->packedMode()) {
    region = regionSparse2->denseVector();
  } else { // packed mode
    for (int j = 0; j < numberNonZero; j++) {
      region[regionIndex[j]] = region2[j];
      region2[j] = 0.0;
    }
  }

  double *solution = workArea2_;
  ftran(region, solution, save);

  // get nonzeros
  numberNonZero = 0;
  if (!regionSparse2->packedMode()) {
    for (int i = 0; i < numberRows_; i++) {
      const double value = solution[i];
      if (fabs(value) > zeroTolerance_) {
        region[i] = value;
        regionIndex[numberNonZero++] = i;
      } else
        region[i] = 0.0;
    }
  } else { // packed mode
    memset(region, 0, numberRows_ * sizeof(double));
    for (int i = 0; i < numberRows_; i++) {
      const double value = solution[i];
      if (fabs(value) > zeroTolerance_) {
        region2[numberNonZero] = value;
        regionIndex[numberNonZero++] = i;
      }
    }
  }
  regionSparse2->setNumElements(numberNonZero);
  return 0;
}

int CoinSimpFactorization::updateTwoColumnsFT(CoinIndexedVector *regionSparse1,
  CoinIndexedVector *regionSparse2,
  CoinIndexedVector *regionSparse3,
  bool)
{
  assert(numberRows_ == numberColumns_);

  double *region2 = regionSparse2->denseVector();
  int *regionIndex2 = regionSparse2->getIndices();
  int numberNonZero2 = regionSparse2->getNumElements();

  double *vec1 = regionSparse1->denseVector();
  if (!regionSparse2->packedMode()) {
    vec1 = regionSparse2->denseVector();
  } else { // packed mode
    for (int j = 0; j < numberNonZero2; j++) {
      vec1[regionIndex2[j]] = region2[j];
      region2[j] = 0.0;
    }
  }
  //
  double *region3 = regionSparse3->denseVector();
  int *regionIndex3 = regionSparse3->getIndices();
  int numberNonZero3 = regionSparse3->getNumElements();
  double *vec2 = auxVector_;
  if (!regionSparse3->packedMode()) {
    vec2 = regionSparse3->denseVector();
  } else { // packed mode
    memset(vec2, 0, numberRows_ * sizeof(double));
    for (int j = 0; j < numberNonZero3; j++) {
      vec2[regionIndex3[j]] = region3[j];
      region3[j] = 0.0;
    }
  }

  double *solution1 = workArea2_;
  double *solution2 = workArea3_;
  ftran2(vec1, solution1, vec2, solution2);

  // get nonzeros
  numberNonZero2 = 0;
  if (!regionSparse2->packedMode()) {
    double value;
    for (int i = 0; i < numberRows_; i++) {
      value = solution1[i];
      if (fabs(value) > zeroTolerance_) {
        vec1[i] = value;
        regionIndex2[numberNonZero2++] = i;
      } else
        vec1[i] = 0.0;
    }
  } else { // packed mode
    double value;
    for (int i = 0; i < numberRows_; i++) {
      vec1[i] = 0.0;
      value = solution1[i];
      if (fabs(value) > zeroTolerance_) {
        region2[numberNonZero2] = value;
        regionIndex2[numberNonZero2++] = i;
      }
    }
  }
  regionSparse2->setNumElements(numberNonZero2);
  //
  numberNonZero3 = 0;
  if (!regionSparse3->packedMode()) {
    double value;
    for (int i = 0; i < numberRows_; i++) {
      value = solution2[i];
      if (fabs(value) > zeroTolerance_) {
        vec2[i] = value;
        regionIndex3[numberNonZero3++] = i;
      } else
        vec2[i] = 0.0;
    }
  } else { // packed mode
    double value;
    for (int i = 0; i < numberRows_; i++) {
      value = solution2[i];
      if (fabs(value) > zeroTolerance_) {
        region3[numberNonZero3] = value;
        regionIndex3[numberNonZero3++] = i;
      }
    }
  }
  regionSparse3->setNumElements(numberNonZero3);
  return 0;
}

int CoinSimpFactorization::updateColumnTranspose(CoinIndexedVector *regionSparse,
  CoinIndexedVector *regionSparse2) const
{
  upColumnTranspose(regionSparse, regionSparse2);
  return 0;
}

int CoinSimpFactorization::upColumnTranspose(CoinIndexedVector *regionSparse,
  CoinIndexedVector *regionSparse2) const
{
  assert(numberRows_ == numberColumns_);
  double *region2 = regionSparse2->denseVector();
  int *regionIndex = regionSparse2->getIndices();
  int numberNonZero = regionSparse2->getNumElements();
  double *region = regionSparse->denseVector();
  if (!regionSparse2->packedMode()) {
    region = regionSparse2->denseVector();
  } else { // packed
    for (int j = 0; j < numberNonZero; j++) {
      region[regionIndex[j]] = region2[j];
      region2[j] = 0.0;
    }
  }
  double *solution = workArea2_;
  btran(region, solution);
  // get nonzeros
  numberNonZero = 0;
  if (!regionSparse2->packedMode()) {
    double value;
    for (int i = 0; i < numberRows_; i++) {
      value = solution[i];
      if (fabs(value) > zeroTolerance_) {
        region[i] = value;
        regionIndex[numberNonZero++] = i;
      } else
        region[i] = 0.0;
    }
  } else { // packed mode
    memset(region, 0, numberRows_ * sizeof(double));
    for (int i = 0; i < numberRows_; i++) {
      const double value = solution[i];
      if (fabs(value) > zeroTolerance_) {
        region2[numberNonZero] = value;
        regionIndex[numberNonZero++] = i;
      }
    }
  }
  regionSparse2->setNumElements(numberNonZero);
  return 0;
}

int CoinSimpFactorization::mainLoopFactor(FactorPointers &pointers)
{
  numberGoodU_ = 0;
  numberSlacks_ = 0;
  bool ifSlack = true;
  for (int i = 0; i < numberColumns_; ++i) {
    int r, s;
    //s=i;
    if (findPivot(pointers, r, s, ifSlack)) {
      return -1;
    }
    if (ifSlack)
      ++numberSlacks_;
    const int rowPos = rowPosition_[r];
    const int colPos = colPosition_[s];
    assert(i <= rowPos && rowPos < numberRows_);
    assert(i <= colPos && colPos < numberColumns_);
    // permute columns
    int j = colOfU_[i];
    colOfU_[i] = colOfU_[colPos];
    colOfU_[colPos] = j;
    colPosition_[colOfU_[i]] = i;
    colPosition_[colOfU_[colPos]] = colPos;
    // permute rows
    j = rowOfU_[i];
    rowOfU_[i] = rowOfU_[rowPos];
    rowOfU_[rowPos] = j;
    rowPosition_[rowOfU_[i]] = i;
    rowPosition_[rowOfU_[rowPos]] = rowPos;
    GaussEliminate(pointers, r, s);
    //if ( maxU_ > maxGrowth_ * maxA_  ){
    //  return -3;
    //}
    ++numberGoodU_;
  }
  return 0;
}
/// find a pivot in the active part of U
int CoinSimpFactorization::findPivot(FactorPointers &pointers, int &r,
  int &s, bool &ifSlack)
{
  int *firstRowKnonzeros = pointers.firstRowKnonzeros;
  int *nextRow = pointers.nextRow;
  int *firstColKnonzeros = pointers.firstColKnonzeros;
  int *prevColumn = pointers.prevColumn;
  int *nextColumn = pointers.nextColumn;
  r = s = -1;
  int numCandidates = 0;
  double bestMarkowitzCount = COIN_DBL_MAX;
  // if there is a column with one element choose it as pivot
  int column = firstColKnonzeros[1];
  if (column != -1) {
    assert(UcolLengths_[column] == 1);
    r = UcolInd_[UcolStarts_[column]];
    s = column;
    if (!colSlack_[column])
      ifSlack = false;
    return 0;
  }
  // from now on no more slacks
  ifSlack = false;
  // if there is a  row with one element choose it
  int row = firstRowKnonzeros[1];
  if (row != -1) {
    assert(UrowLengths_[row] == 1);
    s = UrowInd_[UrowStarts_[row]];
    r = row;
    return 0;
  }
  // consider other rows and columns
  for (int length = 2; length <= numberRows_; ++length) {
    int nextCol = -1;
    for (column = firstColKnonzeros[length]; column != -1; column = nextCol) {
      nextCol = nextColumn[column];
      int minRow, minRowLength;
      int rc = findShortRow(column, length, minRow, minRowLength, pointers);
      if (rc == 0) {
        r = minRow;
        s = column;
        return 0;
      }
      if (minRow != -1) {
        ++numCandidates;
        double MarkowitzCount = static_cast< double >(minRowLength - 1) * (length - 1);
        if (MarkowitzCount < bestMarkowitzCount) {
          r = minRow;
          s = column;
          bestMarkowitzCount = MarkowitzCount;
        }
        if (numCandidates == pivotCandLimit_)
          return 0;
      } else {
        if (doSuhlHeuristic_) {
          // this column did not give a candidate, it will be
          // removed until it becomes a singleton
          removeColumnFromActSet(column, pointers);
          prevColumn[column] = nextColumn[column] = column;
        }
      }
    } // end for ( column= ....
    // now rows
    for (row = firstRowKnonzeros[length]; row != -1; row = nextRow[row]) {
      int minCol, minColLength;
      int rc = findShortColumn(row, length, minCol, minColLength, pointers);
      if (rc == 0) {
        r = row;
        s = minCol;
        return 0;
      }
      if (minCol != -1) {
        ++numCandidates;
        double MarkowitzCount = static_cast< double >(minColLength - 1) * (length - 1);
        if (MarkowitzCount < bestMarkowitzCount) {
          r = row;
          s = minCol;
          bestMarkowitzCount = MarkowitzCount;
        }
        if (numCandidates == pivotCandLimit_)
          return 0;
      }
      //else abort();
    } // end for ( row= ...
  } // end for ( int length= ...
  if (r == -1 || s == -1)
    return 1;
  else
    return 0;
}
//
int CoinSimpFactorization::findPivotShCol(FactorPointers &pointers, int &r, int &s)
{
  int *firstColKnonzeros = pointers.firstColKnonzeros;
  r = s = -1;
  // if there is a column with one element choose it as pivot
  int column = firstColKnonzeros[1];
  if (column != -1) {
    assert(UcolLengths_[column] == 1);
    r = UcolInd_[UcolStarts_[column]];
    s = column;
    return 0;
  }
  // consider other columns
  for (int length = 2; length <= numberRows_; ++length) {
    column = firstColKnonzeros[length];
    if (column != -1)
      break;
  }
  if (column == -1)
    return 1;
  // find largest element
  const int colBeg = UcolStarts_[column];
  const int colEnd = colBeg + UcolLengths_[column];
  double largest = 0.0;
  int rowLargest = -1;
  for (int j = colBeg; j < colEnd; ++j) {
    const int row = UcolInd_[j];
    const int columnIndx = findInRow(row, column);
    assert(columnIndx != -1);
    double coeff = fabs(Urows_[columnIndx]);
    if (coeff < largest)
      continue;
    largest = coeff;
    rowLargest = row;
  }
  assert(rowLargest != -1);
  s = column;
  r = rowLargest;
  return 0;
}

int CoinSimpFactorization::findPivotSimp(FactorPointers &, int &r, int &s)
{
  r = -1;
  int column = s;
  const int colBeg = UcolStarts_[column];
  const int colEnd = colBeg + UcolLengths_[column];
  double largest = 0.0;
  int rowLargest = -1;
  for (int j = colBeg; j < colEnd; ++j) {
    const int row = UcolInd_[j];
    const int columnIndx = findInRow(row, column);
    assert(columnIndx != -1);
    double coeff = fabs(Urows_[columnIndx]);
    if (coeff < largest)
      continue;
    largest = coeff;
    rowLargest = row;
  }
  if (rowLargest != -1) {
    r = rowLargest;
    return 0;
  } else
    return 1;
}

int CoinSimpFactorization::findShortRow(const int column,
  const int length,
  int &minRow,
  int &minRowLength,
  FactorPointers &pointers)
{
  const int colBeg = UcolStarts_[column];
  const int colEnd = colBeg + UcolLengths_[column];
  minRow = -1;
  minRowLength = COIN_INT_MAX;
  for (int j = colBeg; j < colEnd; ++j) {
    int row = UcolInd_[j];
    if (UrowLengths_[row] >= minRowLength)
      continue;
    double largestInRow = findMaxInRrow(row, pointers);
    // find column in row
    int columnIndx = findInRow(row, column);
    assert(columnIndx != -1);
    double coeff = Urows_[columnIndx];
    if (fabs(coeff) < pivotTolerance_ * largestInRow)
      continue;
    minRow = row;
    minRowLength = UrowLengths_[row];
    if (UrowLengths_[row] <= length)
      return 0;
  }
  return 1;
}
int CoinSimpFactorization::findShortColumn(const int row,
  const int length,
  int &minCol,
  int &minColLength,
  FactorPointers &pointers)
{
  const int rowBeg = UrowStarts_[row];
  const int rowEnd = rowBeg + UrowLengths_[row];
  minCol = -1;
  minColLength = COIN_INT_MAX;
  double largestInRow = findMaxInRrow(row, pointers);
  for (int i = rowBeg; i < rowEnd; ++i) {
    int column = UrowInd_[i];
    if (UcolLengths_[column] >= minColLength)
      continue;
    double coeff = Urows_[i];
    if (fabs(coeff) < pivotTolerance_ * largestInRow)
      continue;
    minCol = column;
    minColLength = UcolLengths_[column];
    if (minColLength <= length)
      return 0;
  }
  return 1;
}
// Gaussian elimination
void CoinSimpFactorization::GaussEliminate(FactorPointers &pointers, int &pivotRow, int &pivotCol)
{
  assert(pivotRow >= 0 && pivotRow < numberRows_);
  assert(pivotCol >= 0 && pivotCol < numberRows_);
  int *firstColKnonzeros = pointers.firstColKnonzeros;
  int *prevColumn = pointers.prevColumn;
  int *nextColumn = pointers.nextColumn;
  int *colLabels = vecLabels_;
  double *denseRow = denseVector_;
  removeRowFromActSet(pivotRow, pointers);
  removeColumnFromActSet(pivotCol, pointers);
  // find column s
  int indxColS = findInRow(pivotRow, pivotCol);
  assert(indxColS >= 0);
  // store the inverse of the pivot and remove it from row
  double invPivot = 1.0 / Urows_[indxColS];
  invOfPivots_[pivotRow] = invPivot;
  int rowBeg = UrowStarts_[pivotRow];
  int rowEnd = rowBeg + UrowLengths_[pivotRow];
  Urows_[indxColS] = Urows_[rowEnd - 1];
  UrowInd_[indxColS] = UrowInd_[rowEnd - 1];
  --UrowLengths_[pivotRow];
  --rowEnd;
  // now remove pivot from column
  int indxRowR = findInColumn(pivotCol, pivotRow);
  assert(indxRowR >= 0);
  const int pivColEnd = UcolStarts_[pivotCol] + UcolLengths_[pivotCol];
  UcolInd_[indxRowR] = UcolInd_[pivColEnd - 1];
  --UcolLengths_[pivotCol];
  // go through pivot row
  for (int i = rowBeg; i < rowEnd; ++i) {
    int column = UrowInd_[i];
    colLabels[column] = 1;
    denseRow[column] = Urows_[i];
    // remove this column from bucket because it will change
    removeColumnFromActSet(column, pointers);
    // remove element (pivotRow, column) from column
    int indxRow = findInColumn(column, pivotRow);
    assert(indxRow >= 0);
    const int colEnd = UcolStarts_[column] + UcolLengths_[column];
    UcolInd_[indxRow] = UcolInd_[colEnd - 1];
    --UcolLengths_[column];
  }
  //
  pivoting(pivotRow, pivotCol, invPivot, pointers);
  //
  rowBeg = UrowStarts_[pivotRow];
  rowEnd = rowBeg + UrowLengths_[pivotRow];
  for (int i = rowBeg; i < rowEnd; ++i) {
    int column = UrowInd_[i];
    // clean back these two arrays
    colLabels[column] = 0;
    denseRow[column] = 0.0;
    // column goes into a bucket, if Suhl' heuristic had removed it, it
    // can go back only if it is a singleton
    if (UcolLengths_[column] != 1 || prevColumn[column] != column || nextColumn[column] != column) {
      prevColumn[column] = -1;
      nextColumn[column] = firstColKnonzeros[UcolLengths_[column]];
      if (nextColumn[column] != -1)
        prevColumn[nextColumn[column]] = column;
      firstColKnonzeros[UcolLengths_[column]] = column;
    }
  }
}
void CoinSimpFactorization::pivoting(const int pivotRow,
  const int pivotColumn,
  const double invPivot,
  FactorPointers &pointers)
{
  // initialize the new column of L
  LcolStarts_[pivotRow] = LcolSize_;
  // go trough pivot column
  const int colBeg = UcolStarts_[pivotColumn];
  int colEnd = colBeg + UcolLengths_[pivotColumn];
  for (int i = colBeg; i < colEnd; ++i) {
    int row = UcolInd_[i];
    // remove row from bucket because it will change
    removeRowFromActSet(row, pointers);
    // find pivot column
    int pivotColInRow = findInRow(row, pivotColumn);
    assert(pivotColInRow >= 0);
    const double multiplier = Urows_[pivotColInRow] * invPivot;
    // remove element (row,pivotColumn) from row
    const int currentRowEnd = UrowStarts_[row] + UrowLengths_[row];
    Urows_[pivotColInRow] = Urows_[currentRowEnd - 1];
    UrowInd_[pivotColInRow] = UrowInd_[currentRowEnd - 1];
    --UrowLengths_[row];
    int newNonZeros = UrowLengths_[pivotRow];
    updateCurrentRow(pivotRow, row, multiplier, pointers,
      newNonZeros);
    // store multiplier
    if (LcolSize_ == LcolCap_)
      increaseLsize();
    Lcolumns_[LcolSize_] = multiplier;
    LcolInd_[LcolSize_++] = row;
    ++LcolLengths_[pivotRow];
  }
  // remove elements of pivot column
  UcolLengths_[pivotColumn] = 0;
  // remove pivot column from Ucol_
  if (prevColInU_[pivotColumn] == -1)
    firstColInU_ = nextColInU_[pivotColumn];
  else {
    nextColInU_[prevColInU_[pivotColumn]] = nextColInU_[pivotColumn];
#ifdef COIN_SIMP_CAPACITY
    UcolCapacities_[prevColInU_[pivotColumn]] += UcolCapacities_[pivotColumn];
    UcolCapacities_[pivotColumn] = 0;
#endif
  }
  if (nextColInU_[pivotColumn] == -1)
    lastColInU_ = prevColInU_[pivotColumn];
  else
    prevColInU_[nextColInU_[pivotColumn]] = prevColInU_[pivotColumn];
}

void CoinSimpFactorization::updateCurrentRow(const int pivotRow,
  const int row,
  const double multiplier,
  FactorPointers &pointers,
  int &newNonZeros)
{
  double *rowMax = pointers.rowMax;
  int *firstRowKnonzeros = pointers.firstRowKnonzeros;
  int *prevRow = pointers.prevRow;
  int *nextRow = pointers.nextRow;
  int *colLabels = vecLabels_;
  double *denseRow = denseVector_;
  const int rowBeg = UrowStarts_[row];
  int rowEnd = rowBeg + UrowLengths_[row];
  // treat old nonzeros
  for (int i = rowBeg; i < rowEnd; ++i) {
    const int column = UrowInd_[i];
    if (colLabels[column]) {
      Urows_[i] -= multiplier * denseRow[column];
      const double absNewCoeff = fabs(Urows_[i]);
      colLabels[column] = 0;
      --newNonZeros;
      if (absNewCoeff < zeroTolerance_) {
        // remove it from row
        UrowInd_[i] = UrowInd_[rowEnd - 1];
        Urows_[i] = Urows_[rowEnd - 1];
        --UrowLengths_[row];
        --i;
        --rowEnd;
        // remove it from column
        int indxRow = findInColumn(column, row);
        assert(indxRow >= 0);
        const int colEnd = UcolStarts_[column] + UcolLengths_[column];
        UcolInd_[indxRow] = UcolInd_[colEnd - 1];
        --UcolLengths_[column];
      } else {
        if (maxU_ < absNewCoeff)
          maxU_ = absNewCoeff;
      }
    }
  }
  // now add the new nonzeros to the row
#ifdef COIN_SIMP_CAPACITY
  if (UrowLengths_[row] + newNonZeros > UrowCapacities_[row])
    increaseRowSize(row, UrowLengths_[row] + newNonZeros);
#endif
  const int pivotRowBeg = UrowStarts_[pivotRow];
  const int pivotRowEnd = pivotRowBeg + UrowLengths_[pivotRow];
  int numNew = 0;
  int *newCols = pointers.newCols;
  for (int i = pivotRowBeg; i < pivotRowEnd; ++i) {
    const int column = UrowInd_[i];
    if (colLabels[column]) {
      const double value = -multiplier * denseRow[column];
      const double absNewCoeff = fabs(value);
      if (absNewCoeff >= zeroTolerance_) {
        const int newInd = UrowStarts_[row] + UrowLengths_[row];
        Urows_[newInd] = value;
        UrowInd_[newInd] = column;
        ++UrowLengths_[row];
        newCols[numNew++] = column;
        if (maxU_ < absNewCoeff)
          maxU_ = absNewCoeff;
      }
    } else
      colLabels[column] = 1;
  }
  // add the new nonzeros to the columns
  for (int i = 0; i < numNew; ++i) {
    const int column = newCols[i];
#ifdef COIN_SIMP_CAPACITY
    if (UcolLengths_[column] + 1 > UcolCapacities_[column]) {
      increaseColSize(column, UcolLengths_[column] + 1, false);
    }
#endif
    const int newInd = UcolStarts_[column] + UcolLengths_[column];
    UcolInd_[newInd] = row;
    ++UcolLengths_[column];
  }
  // the row goes to a new bucket
  prevRow[row] = -1;
  nextRow[row] = firstRowKnonzeros[UrowLengths_[row]];
  if (nextRow[row] != -1)
    prevRow[nextRow[row]] = row;
  firstRowKnonzeros[UrowLengths_[row]] = row;
  //
  rowMax[row] = -1.0;
}
#ifdef COIN_SIMP_CAPACITY

void CoinSimpFactorization::increaseRowSize(const int row,
  const int newSize)
{
  assert(newSize > UrowCapacities_[row]);
  const int newNumElements = newSize + minIncrease_;
  if (UrowMaxCap_ < UrowEnd_ + newNumElements) {
    enlargeUrow(UrowEnd_ + newNumElements - UrowMaxCap_);
  }
  int currentCapacity = UrowCapacities_[row];
  memcpy(&Urows_[UrowEnd_], &Urows_[UrowStarts_[row]],
    UrowLengths_[row] * sizeof(double));
  memcpy(&UrowInd_[UrowEnd_], &UrowInd_[UrowStarts_[row]],
    UrowLengths_[row] * sizeof(int));
  UrowStarts_[row] = UrowEnd_;
  UrowCapacities_[row] = newNumElements;
  UrowEnd_ += UrowCapacities_[row];
  if (firstRowInU_ == lastRowInU_)
    return; // only one element
  // remove row from list
  if (prevRowInU_[row] == -1)
    firstRowInU_ = nextRowInU_[row];
  else {
    nextRowInU_[prevRowInU_[row]] = nextRowInU_[row];
    UrowCapacities_[prevRowInU_[row]] += currentCapacity;
  }
  if (nextRowInU_[row] == -1)
    lastRowInU_ = prevRowInU_[row];
  else
    prevRowInU_[nextRowInU_[row]] = prevRowInU_[row];
  // add row at the end of list
  nextRowInU_[lastRowInU_] = row;
  nextRowInU_[row] = -1;
  prevRowInU_[row] = lastRowInU_;
  lastRowInU_ = row;
}
#endif
#ifdef COIN_SIMP_CAPACITY
void CoinSimpFactorization::increaseColSize(const int column,
  const int newSize,
  const bool ifElements)
{
  assert(newSize > UcolCapacities_[column]);
  const int newNumElements = newSize + minIncrease_;
  if (UcolMaxCap_ < UcolEnd_ + newNumElements) {
    enlargeUcol(UcolEnd_ + newNumElements - UcolMaxCap_,
      ifElements);
  }
  int currentCapacity = UcolCapacities_[column];
  memcpy(&UcolInd_[UcolEnd_], &UcolInd_[UcolStarts_[column]],
    UcolLengths_[column] * sizeof(int));
  if (ifElements) {
    memcpy(&Ucolumns_[UcolEnd_], &Ucolumns_[UcolStarts_[column]],
      UcolLengths_[column] * sizeof(double));
  }
  UcolStarts_[column] = UcolEnd_;
  UcolCapacities_[column] = newNumElements;
  UcolEnd_ += UcolCapacities_[column];
  if (firstColInU_ == lastColInU_)
    return; // only one column
  // remove from list
  if (prevColInU_[column] == -1)
    firstColInU_ = nextColInU_[column];
  else {
    nextColInU_[prevColInU_[column]] = nextColInU_[column];
    UcolCapacities_[prevColInU_[column]] += currentCapacity;
  }
  if (nextColInU_[column] == -1)
    lastColInU_ = prevColInU_[column];
  else
    prevColInU_[nextColInU_[column]] = prevColInU_[column];
  // add column at the end
  nextColInU_[lastColInU_] = column;
  nextColInU_[column] = -1;
  prevColInU_[column] = lastColInU_;
  lastColInU_ = column;
}
#endif
double CoinSimpFactorization::findMaxInRrow(const int row,
  FactorPointers &pointers)
{
  double *rowMax = pointers.rowMax;
  double largest = rowMax[row];
  if (largest >= 0.0)
    return largest;
  const int rowBeg = UrowStarts_[row];
  const int rowEnd = rowBeg + UrowLengths_[row];
  for (int i = rowBeg; i < rowEnd; ++i) {
    const double absValue = fabs(Urows_[i]);
    if (absValue > largest)
      largest = absValue;
  }
  rowMax[row] = largest;
  return largest;
}
void CoinSimpFactorization::increaseLsize()
{
  int newcap = LcolCap_ + minIncrease_;

  double *aux = new double[newcap];
  memcpy(aux, Lcolumns_, LcolCap_ * sizeof(double));
  delete[] Lcolumns_;
  Lcolumns_ = aux;

  int *iaux = new int[newcap];
  memcpy(iaux, LcolInd_, LcolCap_ * sizeof(int));
  delete[] LcolInd_;
  LcolInd_ = iaux;

  LcolCap_ = newcap;
}
void CoinSimpFactorization::enlargeUcol(const int numNewElements, const bool ifElements)
{
  int *iaux = new int[UcolMaxCap_ + numNewElements];
  memcpy(iaux, UcolInd_, UcolMaxCap_ * sizeof(int));
  delete[] UcolInd_;
  UcolInd_ = iaux;

  if (ifElements) {
    double *aux = new double[UcolMaxCap_ + numNewElements];
    memcpy(aux, Ucolumns_, UcolMaxCap_ * sizeof(double));
    delete[] Ucolumns_;
    Ucolumns_ = aux;
  }

  UcolMaxCap_ += numNewElements;
}
void CoinSimpFactorization::enlargeUrow(const int numNewElements)
{
  int *iaux = new int[UrowMaxCap_ + numNewElements];
  memcpy(iaux, UrowInd_, UrowMaxCap_ * sizeof(int));
  delete[] UrowInd_;
  UrowInd_ = iaux;

  double *aux = new double[UrowMaxCap_ + numNewElements];
  memcpy(aux, Urows_, UrowMaxCap_ * sizeof(double));
  delete[] Urows_;
  Urows_ = aux;

  UrowMaxCap_ += numNewElements;
}

void CoinSimpFactorization::copyUbyColumns()
{
  memset(UcolLengths_, 0, numberColumns_ * sizeof(int));
  for (int column = 0; column < numberColumns_; ++column) {
    prevColInU_[column] = column - 1;
    nextColInU_[column] = column + 1;
  }
  nextColInU_[numberColumns_ - 1] = -1;
  firstColInU_ = 0;
  lastColInU_ = numberColumns_ - 1;
  //int nonZeros=0;
  //for ( int row=0; row<numberRows_; ++row ){
  //const int rowBeg=UrowStarts_[row];
  //const int rowEnd=rowBeg+UrowLengths_[row];
  //for ( int j=rowBeg; j<rowEnd; ++j )
  //   ++UcolCapacities_[UrowInd_[j]];
  //	nonZeros+=UrowLengths_[row];
  //   }
  //
  //memset(UcolCapacities_, numberRows_, numberColumns_ * sizeof(int));
#ifdef COIN_SIMP_CAPACITY
  for (int i = 0; i < numberColumns_; ++i)
    UcolCapacities_[i] = numberRows_;
#endif
  int k = 0;
  for (int column = 0; column < numberColumns_; ++column) {
    UcolStarts_[column] = k;
    //UcolCapacities_[column]+=minIncrease_;
    //k+=UcolCapacities_[column];
    k += numberRows_;
  }
  UcolEnd_ = k;
  // go through the rows and fill the columns, assume UcolLengths_[]=0
  for (int row = 0; row < numberRows_; ++row) {
    const int rowBeg = UrowStarts_[row];
    int rowEnd = rowBeg + UrowLengths_[row];
    for (int j = rowBeg; j < rowEnd; ++j) {
      // remove els close to zero
      while (fabs(Urows_[j]) < zeroTolerance_) {
        --UrowLengths_[row];
        --rowEnd;
        if (j < rowEnd) {
          Urows_[j] = Urows_[rowEnd];
          UrowInd_[j] = UrowInd_[rowEnd];
        } else
          break;
      }
      if (j == rowEnd)
        continue;
      //
      const int column = UrowInd_[j];
      const int indx = UcolStarts_[column] + UcolLengths_[column];
      Ucolumns_[indx] = Urows_[j];
      UcolInd_[indx] = row;
      ++UcolLengths_[column];
    }
  }
}
void CoinSimpFactorization::copyLbyRows()
{
  int nonZeros = 0;
  memset(LrowLengths_, 0, numberRows_ * sizeof(int));
  for (int column = 0; column < numberRows_; ++column) {
    const int colBeg = LcolStarts_[column];
    const int colEnd = colBeg + LcolLengths_[column];
    for (int j = colBeg; j < colEnd; ++j)
      ++LrowLengths_[LcolInd_[j]];
    nonZeros += LcolLengths_[column];
  }
  //
  LrowSize_ = nonZeros;
  int k = 0;
  for (int row = 0; row < numberRows_; ++row) {
    LrowStarts_[row] = k;
    k += LrowLengths_[row];
  }
#ifdef COIN_SIMP_CAPACITY
  //memset(LrowLengths_,0,numberRows_*sizeof(int));
  // fill the rows
  for (int column = 0; column < numberRows_; ++column) {
    const int colBeg = LcolStarts_[column];
    const int colEnd = colBeg + LcolLengths_[column];
    for (int j = colBeg; j < colEnd; ++j) {
      const int row = LcolInd_[j];
      const int indx = LrowStarts_[row]++;
      Lrows_[indx] = Lcolumns_[j];
      LrowInd_[indx] = column;
    }
  }
  // Put back starts
  k = 0;
  for (int row = 0; row < numberRows_; ++row) {
    int next = LrowStarts_[row];
    LrowStarts_[row] = k;
    k = next;
  }
#else
  memset(LrowLengths_, 0, numberRows_ * sizeof(int));
  // fill the rows
  for (int column = 0; column < numberRows_; ++column) {
    const int colBeg = LcolStarts_[column];
    const int colEnd = colBeg + LcolLengths_[column];
    for (int j = colBeg; j < colEnd; ++j) {
      const int row = LcolInd_[j];
      const int indx = LrowStarts_[row] + LrowLengths_[row];
      Lrows_[indx] = Lcolumns_[j];
      LrowInd_[indx] = column;
      ++LrowLengths_[row];
    }
  }
#endif
}

int CoinSimpFactorization::findInRow(const int row,
  const int column)
{
  const int rowBeg = UrowStarts_[row];
  const int rowEnd = rowBeg + UrowLengths_[row];
  int columnIndx = -1;
  for (int i = rowBeg; i < rowEnd; ++i) {
    if (UrowInd_[i] == column) {
      columnIndx = i;
      break;
    }
  }
  return columnIndx;
}
int CoinSimpFactorization::findInColumn(const int column,
  const int row)
{
  int indxRow = -1;
  const int colBeg = UcolStarts_[column];
  const int colEnd = colBeg + UcolLengths_[column];
  for (int i = colBeg; i < colEnd; ++i) {
    if (UcolInd_[i] == row) {
      indxRow = i;
      break;
    }
  }
  return indxRow;
}
void CoinSimpFactorization::removeRowFromActSet(const int row,
  FactorPointers &pointers)
{
  int *firstRowKnonzeros = pointers.firstRowKnonzeros;
  int *prevRow = pointers.prevRow;
  int *nextRow = pointers.nextRow;
  if (prevRow[row] == -1)
    firstRowKnonzeros[UrowLengths_[row]] = nextRow[row];
  else
    nextRow[prevRow[row]] = nextRow[row];
  if (nextRow[row] != -1)
    prevRow[nextRow[row]] = prevRow[row];
}
void CoinSimpFactorization::removeColumnFromActSet(const int column,
  FactorPointers &pointers)
{
  int *firstColKnonzeros = pointers.firstColKnonzeros;
  int *prevColumn = pointers.prevColumn;
  int *nextColumn = pointers.nextColumn;
  if (prevColumn[column] == -1)
    firstColKnonzeros[UcolLengths_[column]] = nextColumn[column];
  else
    nextColumn[prevColumn[column]] = nextColumn[column];
  if (nextColumn[column] != -1)
    prevColumn[nextColumn[column]] = prevColumn[column];
}
void CoinSimpFactorization::allocateSomeArrays()
{
  if (denseVector_)
    delete[] denseVector_;
  denseVector_ = new double[numberRows_];
  memset(denseVector_, 0, numberRows_ * sizeof(double));
  if (workArea2_)
    delete[] workArea2_;
  workArea2_ = new double[numberRows_];
  if (workArea3_)
    delete[] workArea3_;
  workArea3_ = new double[numberRows_];

  if (vecLabels_)
    delete[] vecLabels_;
  vecLabels_ = new int[numberRows_];
  memset(vecLabels_, 0, numberRows_ * sizeof(int));
  if (indVector_)
    delete[] indVector_;
  indVector_ = new int[numberRows_];

  if (auxVector_)
    delete[] auxVector_;
  auxVector_ = new double[numberRows_];
  if (auxInd_)
    delete[] auxInd_;
  auxInd_ = new int[numberRows_];

  if (vecKeep_)
    delete[] vecKeep_;
  vecKeep_ = new double[numberRows_];
  if (indKeep_)
    delete[] indKeep_;
  indKeep_ = new int[numberRows_];

  if (LrowStarts_)
    delete[] LrowStarts_;
  LrowStarts_ = new int[numberRows_];
  if (LrowLengths_)
    delete[] LrowLengths_;
  LrowLengths_ = new int[numberRows_];

  LrowCap_ = (numberRows_ * (numberRows_ - 1)) / 2;
  if (Lrows_)
    delete[] Lrows_;
  Lrows_ = new double[LrowCap_];
  if (LrowInd_)
    delete[] LrowInd_;
  LrowInd_ = new int[LrowCap_];

  if (LcolStarts_)
    delete[] LcolStarts_;
  LcolStarts_ = new int[numberRows_];
  if (LcolLengths_)
    delete[] LcolLengths_;
  LcolLengths_ = new int[numberRows_];
  LcolCap_ = LrowCap_;
  if (Lcolumns_)
    delete[] Lcolumns_;
  Lcolumns_ = new double[LcolCap_];
  if (LcolInd_)
    delete[] LcolInd_;
  LcolInd_ = new int[LcolCap_];

  if (UrowStarts_)
    delete[] UrowStarts_;
  UrowStarts_ = new int[numberRows_];
  if (UrowLengths_)
    delete[] UrowLengths_;
  UrowLengths_ = new int[numberRows_];
#ifdef COIN_SIMP_CAPACITY
  if (UrowCapacities_)
    delete[] UrowCapacities_;
  UrowCapacities_ = new int[numberRows_];
#endif

  minIncrease_ = 10;
  UrowMaxCap_ = numberRows_ * (numberRows_ + minIncrease_);
  if (Urows_)
    delete[] Urows_;
  Urows_ = new double[UrowMaxCap_];
  if (UrowInd_)
    delete[] UrowInd_;
  UrowInd_ = new int[UrowMaxCap_];

  if (prevRowInU_)
    delete[] prevRowInU_;
  prevRowInU_ = new int[numberRows_];
  if (nextRowInU_)
    delete[] nextRowInU_;
  nextRowInU_ = new int[numberRows_];

  if (UcolStarts_)
    delete[] UcolStarts_;
  UcolStarts_ = new int[numberRows_];
  if (UcolLengths_)
    delete[] UcolLengths_;
  UcolLengths_ = new int[numberRows_];
#ifdef COIN_SIMP_CAPACITY
  if (UcolCapacities_)
    delete[] UcolCapacities_;
  UcolCapacities_ = new int[numberRows_];
#endif

  UcolMaxCap_ = UrowMaxCap_;
  if (Ucolumns_)
    delete[] Ucolumns_;
  Ucolumns_ = new double[UcolMaxCap_];
  if (UcolInd_)
    delete[] UcolInd_;
  UcolInd_ = new int[UcolMaxCap_];

  if (prevColInU_)
    delete[] prevColInU_;
  prevColInU_ = new int[numberRows_];
  if (nextColInU_)
    delete[] nextColInU_;
  nextColInU_ = new int[numberRows_];
  if (colSlack_)
    delete[] colSlack_;
  colSlack_ = new int[numberRows_];

  if (invOfPivots_)
    delete[] invOfPivots_;
  invOfPivots_ = new double[numberRows_];

  if (colOfU_)
    delete[] colOfU_;
  colOfU_ = new int[numberRows_];
  if (colPosition_)
    delete[] colPosition_;
  colPosition_ = new int[numberRows_];
  if (rowOfU_)
    delete[] rowOfU_;
  rowOfU_ = new int[numberRows_];
  if (rowPosition_)
    delete[] rowPosition_;
  rowPosition_ = new int[numberRows_];
  if (secRowOfU_)
    delete[] secRowOfU_;
  secRowOfU_ = new int[numberRows_];
  if (secRowPosition_)
    delete[] secRowPosition_;
  secRowPosition_ = new int[numberRows_];

  if (EtaPosition_)
    delete[] EtaPosition_;
  EtaPosition_ = new int[maximumPivots_];
  if (EtaStarts_)
    delete[] EtaStarts_;
  EtaStarts_ = new int[maximumPivots_];
  if (EtaLengths_)
    delete[] EtaLengths_;
  EtaLengths_ = new int[maximumPivots_];
  maxEtaRows_ = maximumPivots_;

  EtaMaxCap_ = maximumPivots_ * minIncrease_;
  if (EtaInd_)
    delete[] EtaInd_;
  EtaInd_ = new int[EtaMaxCap_];
  if (Eta_)
    delete[] Eta_;
  Eta_ = new double[EtaMaxCap_];
}

void CoinSimpFactorization::Lxeqb(double *b) const
{
  double *rhs = b;
  int k, colBeg, *ind, *indEnd;
  double xk, *Lcol;
  // now solve
  for (int j = firstNumberSlacks_; j < numberRows_; ++j) {
    k = rowOfU_[j];
    xk = rhs[k];
    if (xk != 0.0) {
      //if ( fabs(xk)>zeroTolerance_ ) {
      colBeg = LcolStarts_[k];
      ind = LcolInd_ + colBeg;
      indEnd = ind + LcolLengths_[k];
      Lcol = Lcolumns_ + colBeg;
      for (; ind != indEnd; ++ind) {
        rhs[*ind] -= (*Lcol) * xk;
        ++Lcol;
      }
    }
  }
}

void CoinSimpFactorization::Lxeqb2(double *b1, double *b2) const
{
  double *rhs1 = b1;
  double *rhs2 = b2;
  double x1, x2, *Lcol;
  int k, colBeg, *ind, *indEnd, j;
  // now solve
  for (j = firstNumberSlacks_; j < numberRows_; ++j) {
    k = rowOfU_[j];
    x1 = rhs1[k];
    x2 = rhs2[k];
    if (x1 == 0.0) {
      if (x2 == 0.0) {
      } else {
        colBeg = LcolStarts_[k];
        ind = LcolInd_ + colBeg;
        indEnd = ind + LcolLengths_[k];
        Lcol = Lcolumns_ + colBeg;
        for (; ind != indEnd; ++ind) {
#if 0
		    rhs2[ *ind ]-= (*Lcol) * x2;
#else
          double value = rhs2[*ind];
          rhs2[*ind] = value - (*Lcol) * x2;
#endif
          ++Lcol;
        }
      }
    } else {
      if (x2 == 0.0) {
        colBeg = LcolStarts_[k];
        ind = LcolInd_ + colBeg;
        indEnd = ind + LcolLengths_[k];
        Lcol = Lcolumns_ + colBeg;
        for (; ind != indEnd; ++ind) {
#if 0
		    rhs1[ *ind ]-= (*Lcol) * x1;
#else
          double value = rhs1[*ind];
          rhs1[*ind] = value - (*Lcol) * x1;
#endif
          ++Lcol;
        }
      } else {
        colBeg = LcolStarts_[k];
        ind = LcolInd_ + colBeg;
        indEnd = ind + LcolLengths_[k];
        Lcol = Lcolumns_ + colBeg;
        for (; ind != indEnd; ++ind) {
#if 0
		    rhs1[ *ind ]-= (*Lcol) * x1;
		    rhs2[ *ind ]-= (*Lcol) * x2;
#else
          double value1 = rhs1[*ind];
          rhs1[*ind] = value1 - (*Lcol) * x1;
          double value2 = rhs2[*ind];
          rhs2[*ind] = value2 - (*Lcol) * x2;
#endif
          ++Lcol;
        }
      }
    }
  }
}

void CoinSimpFactorization::Uxeqb(double *b, double *sol) const
{
  double *rhs = b;
  int row, column, colBeg, *ind, *indEnd, k;
  double x, *uCol;
  // now solve
  for (k = numberRows_ - 1; k >= numberSlacks_; --k) {
    row = secRowOfU_[k];
    x = rhs[row];
    column = colOfU_[k];
    if (x != 0.0) {
      //if ( fabs(x) > zeroTolerance_ ) {
      x *= invOfPivots_[row];
      colBeg = UcolStarts_[column];
      ind = UcolInd_ + colBeg;
      indEnd = ind + UcolLengths_[column];
      uCol = Ucolumns_ + colBeg;
      for (; ind != indEnd; ++ind) {
#if 0
		rhs[ *ind ]-= (*uCol) * x;
#else
        double value = rhs[*ind];
        rhs[*ind] = value - (*uCol) * x;
#endif
        ++uCol;
      }
      sol[column] = x;
    } else
      sol[column] = 0.0;
  }
  for (k = numberSlacks_ - 1; k >= 0; --k) {
    row = secRowOfU_[k];
    column = colOfU_[k];
    sol[column] = -rhs[row];
  }
}

void CoinSimpFactorization::Uxeqb2(double *b1, double *sol1, double *b2, double *sol2) const
{
  double *rhs1 = b1;
  double *rhs2 = b2;
  int row, column, colBeg, *ind, *indEnd;
  double x1, x2, *uCol;
  // now solve
  for (int k = numberRows_ - 1; k >= numberSlacks_; --k) {
    row = secRowOfU_[k];
    x1 = rhs1[row];
    x2 = rhs2[row];
    column = colOfU_[k];
    if (x1 == 0.0) {
      if (x2 == 0.0) {
        sol1[column] = 0.0;
        sol2[column] = 0.0;
      } else {
        x2 *= invOfPivots_[row];
        colBeg = UcolStarts_[column];
        ind = UcolInd_ + colBeg;
        indEnd = ind + UcolLengths_[column];
        uCol = Ucolumns_ + colBeg;
        for (; ind != indEnd; ++ind) {
#if 0
		    rhs2[ *ind ]-= (*uCol) * x2;
#else
          double value = rhs2[*ind];
          rhs2[*ind] = value - (*uCol) * x2;
#endif
          ++uCol;
        }
        sol1[column] = 0.0;
        sol2[column] = x2;
      }
    } else {
      if (x2 == 0.0) {
        x1 *= invOfPivots_[row];
        colBeg = UcolStarts_[column];
        ind = UcolInd_ + colBeg;
        indEnd = ind + UcolLengths_[column];
        uCol = Ucolumns_ + colBeg;
        for (; ind != indEnd; ++ind) {
#if 0
		    rhs1[ *ind ]-= (*uCol) * x1;
#else
          double value = rhs1[*ind];
          rhs1[*ind] = value - (*uCol) * x1;
#endif
          ++uCol;
        }
        sol1[column] = x1;
        sol2[column] = 0.0;
      } else {
        x1 *= invOfPivots_[row];
        x2 *= invOfPivots_[row];
        colBeg = UcolStarts_[column];
        ind = UcolInd_ + colBeg;
        indEnd = ind + UcolLengths_[column];
        uCol = Ucolumns_ + colBeg;
        for (; ind != indEnd; ++ind) {
#if 0
		    rhs1[ *ind ]-= (*uCol) * x1;
		    rhs2[ *ind ]-= (*uCol) * x2;
#else
          double value1 = rhs1[*ind];
          rhs1[*ind] = value1 - (*uCol) * x1;
          double value2 = rhs2[*ind];
          rhs2[*ind] = value2 - (*uCol) * x2;
#endif
          ++uCol;
        }
        sol1[column] = x1;
        sol2[column] = x2;
      }
    }
  }
  for (int k = numberSlacks_ - 1; k >= 0; --k) {
    row = secRowOfU_[k];
    column = colOfU_[k];
    sol1[column] = -rhs1[row];
    sol2[column] = -rhs2[row];
  }
}

void CoinSimpFactorization::xLeqb(double *b) const
{
  double *rhs = b;
  int k, *ind, *indEnd, j;
  int colBeg;
  double x, *Lcol;
  // find last nonzero
  int last;
  for (last = numberColumns_ - 1; last >= 0; --last) {
    if (rhs[rowOfU_[last]])
      break;
  }
  // this seems to be faster
  if (last >= 0) {
    for (j = last; j >= firstNumberSlacks_; --j) {
      k = rowOfU_[j];
      x = rhs[k];
      colBeg = LcolStarts_[k];
      ind = LcolInd_ + colBeg;
      indEnd = ind + LcolLengths_[k];
      Lcol = Lcolumns_ + colBeg;
      for (; ind != indEnd; ++ind) {
        x -= (*Lcol) * rhs[*ind];
        ++Lcol;
      }
      rhs[k] = x;
    }
  } // if ( last >= 0 ){
}

void CoinSimpFactorization::xUeqb(double *b, double *sol) const
{
  double *rhs = b;
  int row, col, *ind, *indEnd, k;
  double xr;
  // now solve
#if 1
  int rowBeg;
  double *uRow;
  for (k = 0; k < numberSlacks_; ++k) {
    row = secRowOfU_[k];
    col = colOfU_[k];
    xr = rhs[col];
    if (xr != 0.0) {
      //if ( fabs(xr)> zeroTolerance_ ) {
      xr = -xr;
      rowBeg = UrowStarts_[row];
      ind = UrowInd_ + rowBeg;
      indEnd = ind + UrowLengths_[row];
      uRow = Urows_ + rowBeg;
      for (; ind != indEnd; ++ind) {
        rhs[*ind] -= (*uRow) * xr;
        ++uRow;
      }
      sol[row] = xr;
    } else
      sol[row] = 0.0;
  }
  for (k = numberSlacks_; k < numberRows_; ++k) {
    row = secRowOfU_[k];
    col = colOfU_[k];
    xr = rhs[col];
    if (xr != 0.0) {
      //if ( fabs(xr)> zeroTolerance_ ) {
      xr *= invOfPivots_[row];
      rowBeg = UrowStarts_[row];
      ind = UrowInd_ + rowBeg;
      indEnd = ind + UrowLengths_[row];
      uRow = Urows_ + rowBeg;
      for (; ind != indEnd; ++ind) {
        rhs[*ind] -= (*uRow) * xr;
        ++uRow;
      }
      sol[row] = xr;
    } else
      sol[row] = 0.0;
  }
#else
  for (k = 0; k < numberSlacks_; ++k) {
    row = secRowOfU_[k];
    col = colOfU_[k];
    sol[row] = -rhs[col];
  }
  for (k = numberSlacks_; k < numberRows_; ++k) {
    row = secRowOfU_[k];
    col = colOfU_[k];
    xr = rhs[col];
    int colBeg = UcolStarts_[col];
    ind = UcolInd_ + colBeg;
    indEnd = ind + UcolLengths_[col];
    double *uCol = Ucolumns_ + colBeg;
    for (; ind != indEnd; ++ind, ++uCol) {
      int iRow = *ind;
      double value = sol[iRow];
      double elementValue = *uCol;
      xr -= value * elementValue;
    }
    if (xr != 0.0) {
      xr *= invOfPivots_[row];
      sol[row] = xr;
    } else
      sol[row] = 0.0;
  }
#endif
}

int CoinSimpFactorization::LUupdate(int newBasicCol)
{
  //checkU();
  // recover vector kept in ftran
  double *newColumn = vecKeep_;
  int *indNewColumn = indKeep_;
  int sizeNewColumn = keepSize_;

  // remove elements of new column of U
  const int colBeg = UcolStarts_[newBasicCol];
  const int colEnd = colBeg + UcolLengths_[newBasicCol];
  for (int i = colBeg; i < colEnd; ++i) {
    const int row = UcolInd_[i];
    const int colInRow = findInRow(row, newBasicCol);
    assert(colInRow >= 0);
    // remove from row
    const int rowEnd = UrowStarts_[row] + UrowLengths_[row];
    Urows_[colInRow] = Urows_[rowEnd - 1];
    UrowInd_[colInRow] = UrowInd_[rowEnd - 1];
    --UrowLengths_[row];
  }
  UcolLengths_[newBasicCol] = 0;
  // now add new column to U
  int lastRowInU = -1;
  for (int i = 0; i < sizeNewColumn; ++i) {
    //if ( fabs(newColumn[i]) < zeroTolerance_ ) continue;
    const int row = indNewColumn[i];
    // add to row
#ifdef COIN_SIMP_CAPACITY
    if (UrowLengths_[row] + 1 > UrowCapacities_[row])
      increaseRowSize(row, UrowLengths_[row] + 1);
#endif
    const int rowEnd = UrowStarts_[row] + UrowLengths_[row];
    UrowInd_[rowEnd] = newBasicCol;
    Urows_[rowEnd] = newColumn[i];
    ++UrowLengths_[row];
    if (lastRowInU < secRowPosition_[row])
      lastRowInU = secRowPosition_[row];
  }
  // add to Ucolumns
#ifdef COIN_SIMP_CAPACITY
  if (sizeNewColumn > UcolCapacities_[newBasicCol])
    increaseColSize(newBasicCol, sizeNewColumn, true);
#endif
  memcpy(&Ucolumns_[UcolStarts_[newBasicCol]], &newColumn[0],
    sizeNewColumn * sizeof(double));
  memcpy(&UcolInd_[UcolStarts_[newBasicCol]], &indNewColumn[0],
    sizeNewColumn * sizeof(int));
  UcolLengths_[newBasicCol] = sizeNewColumn;

  const int posNewCol = colPosition_[newBasicCol];
  if (lastRowInU < posNewCol) {
    // matrix is singular
    return 1;
  }
  // permutations
  const int rowInU = secRowOfU_[posNewCol];
  const int colInU = colOfU_[posNewCol];
  for (int i = posNewCol; i < lastRowInU; ++i) {
    int indx = secRowOfU_[i + 1];
    secRowOfU_[i] = indx;
    secRowPosition_[indx] = i;
    int jndx = colOfU_[i + 1];
    colOfU_[i] = jndx;
    colPosition_[jndx] = i;
  }
  secRowOfU_[lastRowInU] = rowInU;
  secRowPosition_[rowInU] = lastRowInU;
  colOfU_[lastRowInU] = colInU;
  colPosition_[colInU] = lastRowInU;
  if (posNewCol < numberSlacks_) {
    if (lastRowInU >= numberSlacks_)
      --numberSlacks_;
    else
      numberSlacks_ = lastRowInU;
  }
  // rowInU will be transformed
  // denseVector_ is assumed to be initialized to zero
  const int rowBeg = UrowStarts_[rowInU];
  const int rowEnd = rowBeg + UrowLengths_[rowInU];
  for (int i = rowBeg; i < rowEnd; ++i) {
    const int column = UrowInd_[i];
    denseVector_[column] = Urows_[i];
    // remove element
    const int indxRow = findInColumn(column, rowInU);
    assert(indxRow >= 0);
    const int colEnd = UcolStarts_[column] + UcolLengths_[column];
    UcolInd_[indxRow] = UcolInd_[colEnd - 1];
    Ucolumns_[indxRow] = Ucolumns_[colEnd - 1];
    --UcolLengths_[column];
  }
  UrowLengths_[rowInU] = 0;
  // rowInU is empty
  // increase Eta by (lastRowInU-posNewCol) elements
  newEta(rowInU, lastRowInU - posNewCol);
  assert(!EtaLengths_[lastEtaRow_]);
  int saveSize = EtaSize_;
  ;
  for (int i = posNewCol; i < lastRowInU; ++i) {
    const int row = secRowOfU_[i];
    const int column = colOfU_[i];
    if (denseVector_[column] == 0.0)
      continue;
    const double multiplier = denseVector_[column] * invOfPivots_[row];
    denseVector_[column] = 0.0;
    const int rowBeg = UrowStarts_[row];
    const int rowEnd = rowBeg + UrowLengths_[row];
#if ARRAY
    for (int j = rowBeg; j < rowEnd; ++j) {
      denseVector_[UrowInd_[j]] -= multiplier * Urows_[j];
    }
#else
    int *ind = UrowInd_ + rowBeg;
    int *indEnd = UrowInd_ + rowEnd;
    double *uRow = Urows_ + rowBeg;
    for (; ind != indEnd; ++ind) {
      denseVector_[*ind] -= multiplier * (*uRow);
      ++uRow;
    }
#endif
    // store multiplier
    Eta_[EtaSize_] = multiplier;
    EtaInd_[EtaSize_++] = row;
  }
  if (EtaSize_ != saveSize)
    EtaLengths_[lastEtaRow_] = EtaSize_ - saveSize;
  else
    --lastEtaRow_;
  // inverse of diagonal
  invOfPivots_[rowInU] = 1.0 / denseVector_[colOfU_[lastRowInU]];
  denseVector_[colOfU_[lastRowInU]] = 0.0;
  // now store row
  int newEls = 0;
  for (int i = lastRowInU + 1; i < numberColumns_; ++i) {
    const int column = colOfU_[i];
    const double coeff = denseVector_[column];
    denseVector_[column] = 0.0;
    if (fabs(coeff) < zeroTolerance_)
      continue;
#ifdef COIN_SIMP_CAPACITY
    if (UcolLengths_[column] + 1 > UcolCapacities_[column]) {
      increaseColSize(column, UcolLengths_[column] + 1, true);
    }
#endif
    const int newInd = UcolStarts_[column] + UcolLengths_[column];
    UcolInd_[newInd] = rowInU;
    Ucolumns_[newInd] = coeff;
    ++UcolLengths_[column];
    workArea2_[newEls] = coeff;
    indVector_[newEls++] = column;
  }
#ifdef COIN_SIMP_CAPACITY
  if (UrowCapacities_[rowInU] < newEls)
    increaseRowSize(rowInU, newEls);
#endif
  const int startRow = UrowStarts_[rowInU];
  memcpy(&Urows_[startRow], &workArea2_[0], newEls * sizeof(double));
  memcpy(&UrowInd_[startRow], &indVector_[0], newEls * sizeof(int));
  UrowLengths_[rowInU] = newEls;
  //
  if (fabs(invOfPivots_[rowInU]) > updateTol_)
    return 2;

  return 0;
}

void CoinSimpFactorization::newEta(int row, int numNewElements)
{
  if (lastEtaRow_ == maxEtaRows_ - 1) {
    int *iaux = new int[maxEtaRows_ + minIncrease_];
    memcpy(iaux, EtaPosition_, maxEtaRows_ * sizeof(int));
    delete[] EtaPosition_;
    EtaPosition_ = iaux;

    int *jaux = new int[maxEtaRows_ + minIncrease_];
    memcpy(jaux, EtaStarts_, maxEtaRows_ * sizeof(int));
    delete[] EtaStarts_;
    EtaStarts_ = jaux;

    int *kaux = new int[maxEtaRows_ + minIncrease_];
    memcpy(kaux, EtaLengths_, maxEtaRows_ * sizeof(int));
    delete[] EtaLengths_;
    EtaLengths_ = kaux;

    maxEtaRows_ += minIncrease_;
  }
  if (EtaSize_ + numNewElements > EtaMaxCap_) {
    int number = CoinMax(EtaSize_ + numNewElements - EtaMaxCap_, minIncrease_);

    int *iaux = new int[EtaMaxCap_ + number];
    memcpy(iaux, EtaInd_, EtaSize_ * sizeof(int));
    delete[] EtaInd_;
    EtaInd_ = iaux;

    double *aux = new double[EtaMaxCap_ + number];
    memcpy(aux, Eta_, EtaSize_ * sizeof(double));
    delete[] Eta_;
    Eta_ = aux;

    EtaMaxCap_ += number;
  }
  EtaPosition_[++lastEtaRow_] = row;
  EtaStarts_[lastEtaRow_] = EtaSize_;
  EtaLengths_[lastEtaRow_] = 0;
}
void CoinSimpFactorization::copyRowPermutations()
{
  memcpy(&secRowOfU_[0], &rowOfU_[0],
    numberRows_ * sizeof(int));
  memcpy(&secRowPosition_[0], &rowPosition_[0],
    numberRows_ * sizeof(int));
}

void CoinSimpFactorization::Hxeqb(double *b) const
{
  double *rhs = b;
  int row, rowBeg, *ind, *indEnd;
  double xr, *eta;
  // now solve
  for (int k = 0; k <= lastEtaRow_; ++k) {
    row = EtaPosition_[k];
    rowBeg = EtaStarts_[k];
    xr = 0.0;
    ind = EtaInd_ + rowBeg;
    indEnd = ind + EtaLengths_[k];
    eta = Eta_ + rowBeg;
    for (; ind != indEnd; ++ind) {
      xr += rhs[*ind] * (*eta);
      ++eta;
    }
    rhs[row] -= xr;
  }
}

void CoinSimpFactorization::Hxeqb2(double *b1, double *b2) const
{
  double *rhs1 = b1;
  double *rhs2 = b2;
  int row, rowBeg, *ind, *indEnd;
  double x1, x2, *eta;
  // now solve
  for (int k = 0; k <= lastEtaRow_; ++k) {
    row = EtaPosition_[k];
    rowBeg = EtaStarts_[k];
    x1 = 0.0;
    x2 = 0.0;
    ind = EtaInd_ + rowBeg;
    indEnd = ind + EtaLengths_[k];
    eta = Eta_ + rowBeg;
    for (; ind != indEnd; ++ind) {
      x1 += rhs1[*ind] * (*eta);
      x2 += rhs2[*ind] * (*eta);
      ++eta;
    }
    rhs1[row] -= x1;
    rhs2[row] -= x2;
  }
}

void CoinSimpFactorization::xHeqb(double *b) const
{
  double *rhs = b;
  int row, rowBeg, *ind, *indEnd;
  double xr, *eta;
  // now solve
  for (int k = lastEtaRow_; k >= 0; --k) {
    row = EtaPosition_[k];
    xr = rhs[row];
    if (xr == 0.0)
      continue;
    //if ( fabs(xr) <= zeroTolerance_ ) continue;
    rowBeg = EtaStarts_[k];
    ind = EtaInd_ + rowBeg;
    indEnd = ind + EtaLengths_[k];
    eta = Eta_ + rowBeg;
    for (; ind != indEnd; ++ind) {
      rhs[*ind] -= xr * (*eta);
      ++eta;
    }
  }
}

void CoinSimpFactorization::ftran(double *b, double *sol, bool save) const
{
  Lxeqb(b);
  Hxeqb(b);
  if (save) {
    // keep vector
    keepSize_ = 0;
    for (int i = 0; i < numberRows_; ++i) {
      if (fabs(b[i]) < zeroTolerance_)
        continue;
      vecKeep_[keepSize_] = b[i];
      indKeep_[keepSize_++] = i;
    }
  }
  Uxeqb(b, sol);
}

void CoinSimpFactorization::ftran2(double *b1, double *sol1, double *b2, double *sol2) const
{
  Lxeqb2(b1, b2);
  Hxeqb2(b1, b2);
  // keep vector
  keepSize_ = 0;
  for (int i = 0; i < numberRows_; ++i) {
    if (fabs(b1[i]) < zeroTolerance_)
      continue;
    vecKeep_[keepSize_] = b1[i];
    indKeep_[keepSize_++] = i;
  }
  Uxeqb2(b1, sol1, b2, sol2);
}

void CoinSimpFactorization::btran(double *b, double *sol) const
{
  xUeqb(b, sol);
  xHeqb(sol);
  xLeqb(sol);
}

#endif

/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
*/