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

CoinFactorization1.cpp

/* $Id: CoinFactorization1.cpp 1215 2009-11-05 11:03:04Z forrest $ */
// Copyright (C) 2002, International Business Machines
// Corporation and others.  All Rights Reserved.
#if defined(_MSC_VER)
// Turn off compiler warning about long names
#  pragma warning(disable:4786)
#endif

#include "CoinUtilsConfig.h"

#include <cassert>
#include "CoinFactorization.hpp"
#include "CoinIndexedVector.hpp"
#include "CoinHelperFunctions.hpp"
#include "CoinPackedMatrix.hpp"
#include <stdio.h>
//:class CoinFactorization.  Deals with Factorization and Updates
//  CoinFactorization.  Constructor
00019 CoinFactorization::CoinFactorization (  )
{
  persistenceFlag_=0;
  gutsOfInitialize(7);
}

/// Copy constructor 
00026 CoinFactorization::CoinFactorization ( const CoinFactorization &other)
{
  persistenceFlag_=0;
  gutsOfInitialize(3);
  persistenceFlag_=other.persistenceFlag_;
  gutsOfCopy(other);
}
/// The real work of constructors etc
00034 void CoinFactorization::gutsOfDestructor(int )
{
  delete [] denseArea_;
  delete [] densePermute_;

  elementU_.conditionalDelete();
  startRowU_.conditionalDelete();
  convertRowToColumnU_.conditionalDelete();
  indexRowU_.conditionalDelete();
  indexColumnU_.conditionalDelete();
  startColumnU_.conditionalDelete();
  elementL_.conditionalDelete();
  indexRowL_.conditionalDelete();
  startColumnL_.conditionalDelete();
  startColumnR_.conditionalDelete();
  numberInRow_.conditionalDelete();
  numberInColumn_.conditionalDelete();
  numberInColumnPlus_.conditionalDelete();
  pivotColumn_.conditionalDelete();
  pivotColumnBack_.conditionalDelete();
  firstCount_.conditionalDelete();
  nextCount_.conditionalDelete();
  lastCount_.conditionalDelete();
  permute_.conditionalDelete();
  permuteBack_.conditionalDelete();
  nextColumn_.conditionalDelete();
  lastColumn_.conditionalDelete();
  nextRow_.conditionalDelete();
  lastRow_.conditionalDelete();
  saveColumn_.conditionalDelete();
  markRow_.conditionalDelete();
  pivotRowL_.conditionalDelete();
  pivotRegion_.conditionalDelete();
  elementByRowL_.conditionalDelete();
  startRowL_.conditionalDelete();
  indexColumnL_.conditionalDelete();
  sparse_.conditionalDelete();
  workArea_.conditionalDelete();
  workArea2_.conditionalDelete();
  numberCompressions_ = 0;
  biggerDimension_ = 0;
  numberRows_ = 0;
  numberRowsExtra_ = 0;
  maximumRowsExtra_ = 0;
  numberColumns_ = 0;
  numberColumnsExtra_ = 0;
  maximumColumnsExtra_ = 0;
  numberGoodU_ = 0;
  numberGoodL_ = 0;
  totalElements_ = 0;
  factorElements_ = 0;
  status_ = -1;
  numberSlacks_ = 0;
  numberU_ = 0;
  maximumU_=0;
  lengthU_ = 0;
  lengthAreaU_ = 0;
  numberL_ = 0;
  baseL_ = 0;
  lengthL_ = 0;
  lengthAreaL_ = 0;
  numberR_ = 0;
  lengthR_ = 0;
  lengthAreaR_ = 0;
  denseArea_=NULL;
  densePermute_=NULL;
  // next two offsets but this makes cleaner
  elementR_=NULL;
  indexRowR_=NULL;
  numberDense_=0;
  //persistenceFlag_=0;
  ////denseThreshold_=0;
  
}
// type - 1 bit tolerances etc, 2 rest
00109 void CoinFactorization::gutsOfInitialize(int type)
{
  if ((type&2)!=0) {
    numberCompressions_ = 0;
    biggerDimension_ = 0;
    numberRows_ = 0;
    numberRowsExtra_ = 0;
    maximumRowsExtra_ = 0;
    numberColumns_ = 0;
    numberColumnsExtra_ = 0;
    maximumColumnsExtra_ = 0;
    numberGoodU_ = 0;
    numberGoodL_ = 0;
    totalElements_ = 0;
    factorElements_ = 0;
    status_ = -1;
    numberPivots_ = 0;
    numberSlacks_ = 0;
    numberU_ = 0;
    maximumU_=0;
    lengthU_ = 0;
    lengthAreaU_ = 0;
    numberL_ = 0;
    baseL_ = 0;
    lengthL_ = 0;
    lengthAreaL_ = 0;
    numberR_ = 0;
    lengthR_ = 0;
    lengthAreaR_ = 0;
    elementR_ = NULL;
    indexRowR_ = NULL;
    // always switch off sparse
    sparseThreshold_=0;
    sparseThreshold2_= 0;
    denseArea_ = NULL;
    densePermute_=NULL;
    numberDense_=0;
    if (!persistenceFlag_) {
      workArea_=CoinFactorizationDoubleArrayWithLength();
      workArea2_=CoinUnsignedIntArrayWithLength();
      pivotColumn_=CoinIntArrayWithLength();
    }
  }
  // after 2 because of persistenceFlag_
  if ((type&1)!=0) {
    areaFactor_ = 0.0;
    pivotTolerance_ = 1.0e-1;
    zeroTolerance_ = 1.0e-13;
#ifndef COIN_FAST_CODE
    slackValue_ = -1.0;
#endif
    messageLevel_=0;
    maximumPivots_=200;
    numberTrials_ = 4;
    relaxCheck_=1.0;
#if DENSE_CODE==1
    denseThreshold_=31;
    denseThreshold_=71;
#else
    denseThreshold_=0;
#endif
    biasLU_=2;
    doForrestTomlin_=true;
    persistenceFlag_=0;
  }
  if ((type&4)!=0) {
    // we need to get 1 element arrays for any with length n+1 !!
    startColumnL_.conditionalNew (1);
    startColumnR_.conditionalNew( 1 );
    startRowU_.conditionalNew(1);
    numberInRow_.conditionalNew(1);
    nextRow_.conditionalNew(1);
    lastRow_.conditionalNew( 1 );
    pivotRegion_.conditionalNew(1);
    permuteBack_.conditionalNew(1);
    permute_.conditionalNew(1);
    pivotColumnBack_.conditionalNew(1);
    startColumnU_.conditionalNew( 1 );
    numberInColumn_.conditionalNew(1);
    numberInColumnPlus_.conditionalNew(1);
    pivotColumn_.conditionalNew(1);
    nextColumn_.conditionalNew(1);
    lastColumn_.conditionalNew(1);
    collectStatistics_=false;
    
    // Below are all to collect
    ftranCountInput_=0.0;
    ftranCountAfterL_=0.0;
    ftranCountAfterR_=0.0;
    ftranCountAfterU_=0.0;
    btranCountInput_=0.0;
    btranCountAfterU_=0.0;
    btranCountAfterR_=0.0;
    btranCountAfterL_=0.0;
    
    // We can roll over factorizations
    numberFtranCounts_=0;
    numberBtranCounts_=0;
    
    // While these are averages collected over last 
    ftranAverageAfterL_=0;
    ftranAverageAfterR_=0;
    ftranAverageAfterU_=0;
    btranAverageAfterU_=0;
    btranAverageAfterR_=0;
    btranAverageAfterL_=0; 
#ifdef ZEROFAULT
    startColumnL_.array()[0] = 0;
    startColumnR_.array()[0] = 0;
    startRowU_.array()[0] = 0;
    numberInRow_.array()[0] = 0;
    nextRow_.array()[0] = 0;
    lastRow_.array()[0] = 0;
    pivotRegion_.array()[0] = 0.0;
    permuteBack_.array()[0] = 0;
    permute_.array()[0] = 0;
    pivotColumnBack_.array()[0] = 0;
    startColumnU_.array()[0] = 0;
    numberInColumn_.array()[0] = 0;
    numberInColumnPlus_.array()[0] = 0;
    pivotColumn_.array()[0] = 0;
    nextColumn_.array()[0] = 0;
    lastColumn_.array()[0] = 0;
#endif
  }
}
//Part of LP
00236 int CoinFactorization::factorize (
                         const CoinPackedMatrix & matrix,
                         int rowIsBasic[],
                         int columnIsBasic[],
                         double areaFactor )
{
  // maybe for speed will be better to leave as many regions as possible
  gutsOfDestructor();
  gutsOfInitialize(2);
  // ? is this correct
  //if (biasLU_==2)
  //biasLU_=3;
  if (areaFactor)
    areaFactor_ = areaFactor;
  const int * row = matrix.getIndices();
  const CoinBigIndex * columnStart = matrix.getVectorStarts();
  const int * columnLength = matrix.getVectorLengths(); 
  const double * element = matrix.getElements();
  int numberRows=matrix.getNumRows();
  int numberColumns=matrix.getNumCols();
  int numberBasic = 0;
  CoinBigIndex numberElements=0;
  int numberRowBasic=0;

  // compute how much in basis

  int i;

  for (i=0;i<numberRows;i++) {
    if (rowIsBasic[i]>=0)
      numberRowBasic++;
  }

  numberBasic = numberRowBasic;

  for (i=0;i<numberColumns;i++) {
    if (columnIsBasic[i]>=0) {
      numberBasic++;
      numberElements += columnLength[i];
    }
  }
  if ( numberBasic > numberRows ) {
    return -2; // say too many in basis
  }
  numberElements = 3 * numberBasic + 3 * numberElements + 20000;
  getAreas ( numberRows, numberBasic, numberElements,
           2 * numberElements );
  //fill
  //copy
  numberBasic=0;
  numberElements=0;
  int * indexColumnU = indexColumnU_.array();
  int * indexRowU = indexRowU_.array();
  CoinFactorizationDouble * elementU = elementU_.array();
  for (i=0;i<numberRows;i++) {
    if (rowIsBasic[i]>=0) {
      indexRowU[numberElements]=i;
      indexColumnU[numberElements]=numberBasic;
      elementU[numberElements++]=slackValue_;
      numberBasic++;
    }
  }
  for (i=0;i<numberColumns;i++) {
    if (columnIsBasic[i]>=0) {
      CoinBigIndex j;
      for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
      indexRowU[numberElements]=row[j];
      indexColumnU[numberElements]=numberBasic;
      elementU[numberElements++]=element[j];
      }
      numberBasic++;
    }
  }
  lengthU_ = numberElements;
  maximumU_ = numberElements;

  preProcess ( 0 );
  factor (  );
  numberBasic=0;
  if (status_ == 0) {
    int * permuteBack = permuteBack_.array();
    int * back = pivotColumnBack();
    for (i=0;i<numberRows;i++) {
      if (rowIsBasic[i]>=0) {
      rowIsBasic[i]=permuteBack[back[numberBasic++]];
      }
    }
    for (i=0;i<numberColumns;i++) {
      if (columnIsBasic[i]>=0) {
      columnIsBasic[i]=permuteBack[back[numberBasic++]];
      }
    }
    // Set up permutation vector
    // these arrays start off as copies of permute
    // (and we could use permute_ instead of pivotColumn (not back though))
    CoinMemcpyN ( permute_.array(), numberRows_ , pivotColumn_.array()  );
    CoinMemcpyN ( permuteBack_.array(), numberRows_ , pivotColumnBack()  );
  } else if (status_ == -1) {
    const int * pivotColumn = pivotColumn_.array();
    // mark as basic or non basic
    for (i=0;i<numberRows_;i++) {
      if (rowIsBasic[i]>=0) {
      if (pivotColumn[numberBasic]>=0) 
        rowIsBasic[i]=pivotColumn[numberBasic];
      else 
        rowIsBasic[i]=-1;
      numberBasic++;
      }
    }
    for (i=0;i<numberColumns;i++) {
      if (columnIsBasic[i]>=0) {
      if (pivotColumn[numberBasic]>=0) 
        columnIsBasic[i]=pivotColumn[numberBasic];
       else 
        columnIsBasic[i]=-1;
      numberBasic++;
      }
    }
  }

  return status_;
}
//Given as triplets
00359 int CoinFactorization::factorize (
                       int numberOfRows,
                       int numberOfColumns,
                       CoinBigIndex numberOfElements,
                       CoinBigIndex maximumL,
                       CoinBigIndex maximumU,
                       const int indicesRow[],
                       const int indicesColumn[],
                       const double elements[] ,
                       int permutation[],
                       double areaFactor)
{
  gutsOfDestructor();
  gutsOfInitialize(2);
  if (areaFactor)
    areaFactor_ = areaFactor;
  getAreas ( numberOfRows, numberOfColumns, maximumL, maximumU );
  //copy
  CoinMemcpyN ( indicesRow, numberOfElements, indexRowU_.array() );
  CoinMemcpyN ( indicesColumn, numberOfElements, indexColumnU_.array() );
  int i;
  CoinFactorizationDouble * elementU = elementU_.array();
  for (i=0;i<numberOfElements;i++)
    elementU[i] = elements[i];
  lengthU_ = numberOfElements;
  maximumU_ = numberOfElements;
  preProcess ( 0 );
  factor (  );
  //say which column is pivoting on which row
  if (status_ == 0) {
    int * permuteBack = permuteBack_.array();
    int * back = pivotColumnBack();
    // permute so slacks on own rows etc
    for (i=0;i<numberOfColumns;i++) {
      permutation[i]=permuteBack[back[i]];
    }
    // Set up permutation vector
    // these arrays start off as copies of permute
    // (and we could use permute_ instead of pivotColumn (not back though))
    CoinMemcpyN ( permute_.array(), numberRows_ , pivotColumn_.array()  );
    CoinMemcpyN ( permuteBack_.array(), numberRows_ , pivotColumnBack()  );
  } else if (status_ == -1) {
    const int * pivotColumn = pivotColumn_.array();
    // mark as basic or non basic
    for (i=0;i<numberOfColumns;i++) {
      if (pivotColumn[i]>=0) {
      permutation[i]=pivotColumn[i];
      } else {
      permutation[i]=-1;
      }
    }
  }

  return status_;
}
/* Two part version for flexibility
   This part creates arrays for user to fill.
   maximumL is guessed maximum size of L part of
   final factorization, maximumU of U part.  These are multiplied by
   areaFactor which can be computed by user or internally.  
   returns 0 -okay, -99 memory */
int 
00421 CoinFactorization::factorizePart1 ( int numberOfRows,
                           int ,
                           CoinBigIndex numberOfElements,
                           int * indicesRow[],
                           int * indicesColumn[],
                           CoinFactorizationDouble * elements[],
                           double areaFactor)
{
  // maybe for speed will be better to leave as many regions as possible
  gutsOfDestructor();
  gutsOfInitialize(2);
  if (areaFactor)
    areaFactor_ = areaFactor;
  CoinBigIndex numberElements = 3 * numberOfRows + 3 * numberOfElements + 20000;
  getAreas ( numberOfRows, numberOfRows, numberElements,
           2 * numberElements );
  // need to trap memory for -99 code
  *indicesRow = indexRowU_.array() ;
  *indicesColumn = indexColumnU_.array() ;
  *elements = elementU_.array() ;
  lengthU_ = numberOfElements;
  maximumU_ = numberElements;
  return 0;
}
/* This is part two of factorization
   Arrays belong to factorization and were returned by part 1
   If status okay, permutation has pivot rows.
   If status is singular, then basic variables have +1 and ones thrown out have -COIN_INT_MAX
   to say thrown out.
   returns 0 -okay, -1 singular, -99 memory */
int 
00452 CoinFactorization::factorizePart2 (int permutation[],int exactNumberElements)
{
  lengthU_ = exactNumberElements;
  preProcess ( 0 );
  factor (  );
  //say which column is pivoting on which row
  int i;
  int * permuteBack = permuteBack_.array();
  int * back = pivotColumnBack();
  // permute so slacks on own rows etc
  for (i=0;i<numberColumns_;i++) {
    permutation[i]=permuteBack[back[i]];
  }
  if (status_ == 0) {
    // Set up permutation vector
    // these arrays start off as copies of permute
    // (and we could use permute_ instead of pivotColumn (not back though))
    CoinMemcpyN ( permute_.array(), numberRows_ , pivotColumn_.array()  );
    CoinMemcpyN ( permuteBack_.array(), numberRows_ , pivotColumnBack()  );
  } else if (status_ == -1) {
    const int * pivotColumn = pivotColumn_.array();
    // mark as basic or non basic
    for (i=0;i<numberColumns_;i++) {
      if (pivotColumn[i]>=0) {
      permutation[i]=pivotColumn[i];
      } else {
      permutation[i]=-1;
      }
    }
  }

  return status_;
}

//  ~CoinFactorization.  Destructor
00487 CoinFactorization::~CoinFactorization (  )
{
  gutsOfDestructor();
}

//  show_self.  Debug show object
void
00494 CoinFactorization::show_self (  ) const
{
  int i;

  const int * pivotColumn = pivotColumn_.array();
  for ( i = 0; i < numberRows_; i++ ) {
    std::cout << "r " << i << " " << pivotColumn[i];
    if (pivotColumnBack()) std::cout<< " " << pivotColumnBack()[i];
    std::cout<< " " << permute_.array()[i];
    if (permuteBack_.array()) std::cout<< " " << permuteBack_.array()[i];
    std::cout<< " " << pivotRegion_.array()[i];
    std::cout << std::endl;
  }
  for ( i = 0; i < numberRows_; i++ ) {
    std::cout << "u " << i << " " << numberInColumn_.array()[i] << std::endl;
    int j;
    CoinSort_2(indexRowU_.array()+startColumnU_.array()[i],
             indexRowU_.array()+startColumnU_.array()[i]+numberInColumn_.array()[i],
             elementU_.array()+startColumnU_.array()[i]);
    for ( j = startColumnU_.array()[i]; j < startColumnU_.array()[i] + numberInColumn_.array()[i];
        j++ ) {
      assert (indexRowU_.array()[j]>=0&&indexRowU_.array()[j]<numberRows_);
      assert (elementU_.array()[j]>-1.0e50&&elementU_.array()[j]<1.0e50);
      std::cout << indexRowU_.array()[j] << " " << elementU_.array()[j] << std::endl;
    }
  }
  for ( i = 0; i < numberRows_; i++ ) {
    std::cout << "l " << i << " " << startColumnL_.array()[i + 1] -
      startColumnL_.array()[i] << std::endl;
    CoinSort_2(indexRowL_.array()+startColumnL_.array()[i],
             indexRowL_.array()+startColumnL_.array()[i+1],
             elementL_.array()+startColumnL_.array()[i]);
    int j;
    for ( j = startColumnL_.array()[i]; j < startColumnL_.array()[i + 1]; j++ ) {
      std::cout << indexRowL_.array()[j] << " " << elementL_.array()[j] << std::endl;
    }
  }

}
//  sort so can compare
void
00535 CoinFactorization::sort (  ) const
{
  int i;

  for ( i = 0; i < numberRows_; i++ ) {
    CoinSort_2(indexRowU_.array()+startColumnU_.array()[i],
             indexRowU_.array()+startColumnU_.array()[i]+numberInColumn_.array()[i],
             elementU_.array()+startColumnU_.array()[i]);
  }
  for ( i = 0; i < numberRows_; i++ ) {
    CoinSort_2(indexRowL_.array()+startColumnL_.array()[i],
             indexRowL_.array()+startColumnL_.array()[i+1],
             elementL_.array()+startColumnL_.array()[i]);
  }

}

//  getAreas.  Gets space for a factorization
//called by constructors
void
00555 CoinFactorization::getAreas ( int numberOfRows,
                   int numberOfColumns,
                   CoinBigIndex maximumL,
                   CoinBigIndex maximumU )
{
  
  numberRows_ = numberOfRows;
  numberColumns_ = numberOfColumns;
  maximumRowsExtra_ = numberRows_ + maximumPivots_;
  numberRowsExtra_ = numberRows_;
  maximumColumnsExtra_ = numberColumns_ + maximumPivots_;
  numberColumnsExtra_ = numberColumns_;
  lengthAreaU_ = maximumU;
  lengthAreaL_ = maximumL;
  if ( !areaFactor_ ) {
    areaFactor_ = 1.0;
  }
  if ( areaFactor_ != 1.0 ) {
    if ((messageLevel_&16)!=0) 
      printf("Increasing factorization areas by %g\n",areaFactor_);
    lengthAreaU_ =  static_cast<CoinBigIndex> (areaFactor_*lengthAreaU_);
    lengthAreaL_ =  static_cast<CoinBigIndex> (areaFactor_*lengthAreaL_);
  }
  elementU_.conditionalNew( lengthAreaU_ );
  indexRowU_.conditionalNew( lengthAreaU_ );
  indexColumnU_.conditionalNew( lengthAreaU_ );
  elementL_.conditionalNew( lengthAreaL_ );
  indexRowL_.conditionalNew( lengthAreaL_ );
  if (persistenceFlag_) {
    // But we can use all we have if bigger
    int length;
    length = CoinMin(elementU_.getSize(),indexRowU_.getSize());
    if (length>lengthAreaU_) {
      lengthAreaU_=length;
      assert (indexColumnU_.getSize()==indexRowU_.getSize());
    }
    length = CoinMin(elementL_.getSize(),indexRowL_.getSize());
    if (length>lengthAreaL_) {
      lengthAreaL_=length;
    }
  }
  startColumnL_.conditionalNew( numberRows_ + 1 );
  startColumnL_.array()[0] = 0;
  startRowU_.conditionalNew( maximumRowsExtra_ + 1);
  // make sure this is valid
  startRowU_.array()[maximumRowsExtra_]=0;
  numberInRow_.conditionalNew( maximumRowsExtra_ + 1 );
  markRow_.conditionalNew( numberRows_ );
  pivotRowL_.conditionalNew( numberRows_ + 1 );
  nextRow_.conditionalNew( maximumRowsExtra_ + 1 );
  lastRow_.conditionalNew( maximumRowsExtra_ + 1 );
  permute_.conditionalNew( maximumRowsExtra_ + 1 );
  pivotRegion_.conditionalNew( maximumRowsExtra_ + 1 );
#ifdef ZEROFAULT
  memset(elementU_.array(),'a',lengthAreaU_*sizeof(CoinFactorizationDouble));
  memset(indexRowU_.array(),'b',lengthAreaU_*sizeof(int));
  memset(indexColumnU_.array(),'c',lengthAreaU_*sizeof(int));
  memset(elementL_.array(),'d',lengthAreaL_*sizeof(CoinFactorizationDouble));
  memset(indexRowL_.array(),'e',lengthAreaL_*sizeof(int));
  memset(startColumnL_.array()+1,'f',numberRows_*sizeof(CoinBigIndex));
  memset(startRowU_.array(),'g',maximumRowsExtra_*sizeof(CoinBigIndex));
  memset(numberInRow_.array(),'h',(maximumRowsExtra_+1)*sizeof(int));
  memset(markRow_.array(),'i',numberRows_*sizeof(int));
  memset(pivotRowL_.array(),'j',(numberRows_+1)*sizeof(int));
  memset(nextRow_.array(),'k',(maximumRowsExtra_+1)*sizeof(int));
  memset(lastRow_.array(),'l',(maximumRowsExtra_+1)*sizeof(int));
  memset(permute_.array(),'l',(maximumRowsExtra_+1)*sizeof(int));
  memset(pivotRegion_.array(),'m',(maximumRowsExtra_+1)*sizeof(CoinFactorizationDouble));
#endif
  startColumnU_.conditionalNew( maximumColumnsExtra_ + 1 );
  numberInColumn_.conditionalNew( maximumColumnsExtra_ + 1 );
  numberInColumnPlus_.conditionalNew( maximumColumnsExtra_ + 1 );
  pivotColumn_.conditionalNew( maximumColumnsExtra_ + 1 );
  nextColumn_.conditionalNew( maximumColumnsExtra_ + 1 );
  lastColumn_.conditionalNew( maximumColumnsExtra_ + 1 );
  saveColumn_.conditionalNew( numberColumns_);
#ifdef ZEROFAULT
  memset(startColumnU_.array(),'a',(maximumColumnsExtra_+1)*sizeof(CoinBigIndex));
  memset(numberInColumn_.array(),'b',(maximumColumnsExtra_+1)*sizeof(int));
  memset(numberInColumnPlus_.array(),'c',(maximumColumnsExtra_+1)*sizeof(int));
  memset(pivotColumn_.array(),'d',(maximumColumnsExtra_+1)*sizeof(int));
  memset(nextColumn_.array(),'e',(maximumColumnsExtra_+1)*sizeof(int));
  memset(lastColumn_.array(),'f',(maximumColumnsExtra_+1)*sizeof(int));
#endif
  if ( numberRows_ + numberColumns_ ) {
    if ( numberRows_ > numberColumns_ ) {
      biggerDimension_ = numberRows_;
    } else {
      biggerDimension_ = numberColumns_;
    }
    firstCount_.conditionalNew( CoinMax(biggerDimension_ + 2, maximumRowsExtra_+1) );
    nextCount_.conditionalNew( numberRows_ + numberColumns_ );
    lastCount_.conditionalNew( numberRows_ + numberColumns_ );
#ifdef ZEROFAULT
    memset(firstCount_.array(),'g',(biggerDimension_ + 2 )*sizeof(int));
    memset(nextCount_.array(),'h',(numberRows_+numberColumns_)*sizeof(int));
    memset(lastCount_.array(),'i',(numberRows_+numberColumns_)*sizeof(int));
#endif
  } else {
    firstCount_.conditionalNew( 2 );
    nextCount_.conditionalNew( 0 );
    lastCount_.conditionalNew( 0 );
#ifdef ZEROFAULT
    memset(firstCount_.array(),'g', 2 *sizeof(int));
#endif
    biggerDimension_ = 0;
  }
}

//  preProcess.  PreProcesses raw triplet data
//state is 0 - triplets, 1 - some counts etc , 2 - ..
void
00667 CoinFactorization::preProcess ( int state,
                     int )
{
  int *indexRow = indexRowU_.array();
  int *indexColumn = indexColumnU_.array();
  CoinFactorizationDouble *element = elementU_.array();
  CoinBigIndex numberElements = lengthU_;
  int *numberInRow = numberInRow_.array();
  int *numberInColumn = numberInColumn_.array();
  int *numberInColumnPlus = numberInColumnPlus_.array();
  CoinBigIndex *startRow = startRowU_.array();
  CoinBigIndex *startColumn = startColumnU_.array();
  int numberRows = numberRows_;
  int numberColumns = numberColumns_;

  if (state<4)
    totalElements_ = numberElements;
  //state falls through to next state
  switch ( state ) {
  case 0:               //counts
    {
      CoinZeroN ( numberInRow, numberRows + 1 );
      CoinZeroN ( numberInColumn, maximumColumnsExtra_ + 1 );
      CoinBigIndex i;
      for ( i = 0; i < numberElements; i++ ) {
      int iRow = indexRow[i];
      int iColumn = indexColumn[i];

      numberInRow[iRow]++;
      numberInColumn[iColumn]++;
      }
    }
  case -1:              //sort
  case 1:               //sort
    {
      CoinBigIndex i, k;

      i = 0;
      int iColumn;
      for ( iColumn = 0; iColumn < numberColumns; iColumn++ ) {
      //position after end of Column
      i += numberInColumn[iColumn];
      startColumn[iColumn] = i;
      }
      for ( k = numberElements - 1; k >= 0; k-- ) {
      int iColumn = indexColumn[k];

      if ( iColumn >= 0 ) {
        CoinFactorizationDouble value = element[k];
        int iRow = indexRow[k];

        indexColumn[k] = -1;
        while ( true ) {
          CoinBigIndex iLook = startColumn[iColumn] - 1;

          startColumn[iColumn] = iLook;
          CoinFactorizationDouble valueSave = element[iLook];
          int iColumnSave = indexColumn[iLook];
          int iRowSave = indexRow[iLook];

          element[iLook] = value;
          indexRow[iLook] = iRow;
          indexColumn[iLook] = -1;
          if ( iColumnSave >= 0 ) {
            iColumn = iColumnSave;
            value = valueSave;
            iRow = iRowSave;
          } else {
            break;
          }
        }               /* endwhile */
      }
      }
    }
  case 2:               //move largest in column to beginning
    //and do row part
    {
      CoinBigIndex i, k;

      i = 0;
      int iRow;
      for ( iRow = 0; iRow < numberRows; iRow++ ) {
      startRow[iRow] = i;
      i += numberInRow[iRow];
      }
      CoinZeroN ( numberInRow, numberRows );
      int iColumn;
      for ( iColumn = 0; iColumn < numberColumns; iColumn++ ) {
      int number = numberInColumn[iColumn];

      if ( number ) {
        CoinBigIndex first = startColumn[iColumn];
        CoinBigIndex largest = first;
        int iRowSave = indexRow[first];
        CoinFactorizationDouble valueSave = element[first];
        double valueLargest = fabs ( valueSave );
        int iLook = numberInRow[iRowSave];

        numberInRow[iRowSave] = iLook + 1;
        indexColumn[startRow[iRowSave] + iLook] = iColumn;
        for ( k = first + 1; k < first + number; k++ ) {
          int iRow = indexRow[k];
          int iLook = numberInRow[iRow];

          numberInRow[iRow] = iLook + 1;
          indexColumn[startRow[iRow] + iLook] = iColumn;
          CoinFactorizationDouble value = element[k];
          double valueAbs = fabs ( value );

          if ( valueAbs > valueLargest ) {
            valueLargest = valueAbs;
            largest = k;
          }
        }
        indexRow[first] = indexRow[largest];
        element[first] = element[largest];
        indexRow[largest] = iRowSave;
        element[largest] = valueSave;
      }
      }
    }
  case 3:               //links and initialize pivots
    {
      int *lastRow = lastRow_.array();
      int *nextRow = nextRow_.array();
      int *lastColumn = lastColumn_.array();
      int *nextColumn = nextColumn_.array();

      CoinFillN ( firstCount_.array(), biggerDimension_ + 2, -1 );
      CoinFillN ( pivotColumn_.array(), numberColumns_, -1 );
      CoinZeroN ( numberInColumnPlus, maximumColumnsExtra_ + 1 );
      int iRow;
      for ( iRow = 0; iRow < numberRows; iRow++ ) {
      lastRow[iRow] = iRow - 1;
      nextRow[iRow] = iRow + 1;
      int number = numberInRow[iRow];

      addLink ( iRow, number );
      }
      lastRow[maximumRowsExtra_] = numberRows - 1;
      nextRow[maximumRowsExtra_] = 0;
      lastRow[0] = maximumRowsExtra_;
      nextRow[numberRows - 1] = maximumRowsExtra_;
      startRow[maximumRowsExtra_] = numberElements;
      int iColumn;
      for ( iColumn = 0; iColumn < numberColumns; iColumn++ ) {
      lastColumn[iColumn] = iColumn - 1;
      nextColumn[iColumn] = iColumn + 1;
      int number = numberInColumn[iColumn];

      addLink ( iColumn + numberRows, number );
      }
      lastColumn[maximumColumnsExtra_] = numberColumns - 1;
      nextColumn[maximumColumnsExtra_] = 0;
      lastColumn[0] = maximumColumnsExtra_;
      if (numberColumns)
      nextColumn[numberColumns - 1] = maximumColumnsExtra_;
      startColumn[maximumColumnsExtra_] = numberElements;
    }
    break;
  case 4:               //move largest in column to beginning
    {
      CoinBigIndex i, k;
      CoinFactorizationDouble * pivotRegion = pivotRegion_.array();
      int iColumn;
      int iRow;
      for ( iRow = 0; iRow < numberRows; iRow++ ) {
      if( numberInRow[iRow]>=0) {
        // zero count
        numberInRow[iRow]=0;
      } else {
        // empty
        //numberInRow[iRow]=-1; already that
      }
      }
      //CoinZeroN ( numberInColumnPlus, maximumColumnsExtra_ + 1 );
      for ( iColumn = 0; iColumn < numberColumns; iColumn++ ) {
      int number = numberInColumn[iColumn];

      if ( number ) {
        // use pivotRegion and startRow for remaining elements
        CoinBigIndex first = startColumn[iColumn];
        CoinBigIndex largest = -1;
        
        double valueLargest = -1.0;
        int nOther=0;
        k = first;
        CoinBigIndex end = first+number;
        for (  ; k < end; k++ ) {
          int iRow = indexRow[k];
          assert (iRow<numberRows_);
          CoinFactorizationDouble value = element[k];
          if (numberInRow[iRow]>=0) {
            numberInRow[iRow]++;
            double valueAbs = fabs ( value );
            if ( valueAbs > valueLargest ) {
            valueLargest = valueAbs;
            largest = nOther;
            }
            startRow[nOther]=iRow;
            pivotRegion[nOther++]=value;
          } else {
            indexRow[first] = iRow;
            element[first++] = value;
          }
        }
        numberInColumnPlus[iColumn]=first-startColumn[iColumn];
        startColumn[iColumn]=first;
        //largest
        if (largest>=0) {
          indexRow[first] = startRow[largest];
          element[first++] = pivotRegion[largest];
        }
        for (k=0;k<nOther;k++) {
          if (k!=largest) {
            indexRow[first] = startRow[k];
            element[first++] = pivotRegion[k];
          }
        }
        numberInColumn[iColumn]=first-startColumn[iColumn];
      }
      }
      //and do row part
      i = 0;
      for ( iRow = 0; iRow < numberRows; iRow++ ) {
      startRow[iRow] = i;
      int n=numberInRow[iRow];
      if (n>0) {
        numberInRow[iRow]=0;
        i += n;
      }
      }
      for ( iColumn = 0; iColumn < numberColumns; iColumn++ ) {
      int number = numberInColumn[iColumn];

      if ( number ) {
        CoinBigIndex first = startColumn[iColumn];
        for ( k = first; k < first + number; k++ ) {
          int iRow = indexRow[k];
          int iLook = numberInRow[iRow];

          numberInRow[iRow] = iLook + 1;
          indexColumn[startRow[iRow] + iLook] = iColumn;
        }
      }
      }
    }
    // modified 3
    {
      //set markRow so no rows updated
      //CoinFillN ( markRow_.array(), numberRows_, -1 );
      int *lastColumn = lastColumn_.array();
      int *nextColumn = nextColumn_.array();
      CoinFactorizationDouble * pivotRegion = pivotRegion_.array();

      int iRow;
      int numberGood=0;
      startColumnL_.array()[0] = 0; //for luck
      for ( iRow = 0; iRow < numberRows; iRow++ ) {
      int number = numberInRow[iRow];
      if (number<0) {
        numberInRow[iRow]=0;
        pivotRegion[numberGood++]=slackValue_;
      }
      }
      int iColumn;
      for ( iColumn = 0 ; iColumn < numberColumns_; iColumn++ ) {
      lastColumn[iColumn] = iColumn - 1;
      nextColumn[iColumn] = iColumn + 1;
      int number = numberInColumn[iColumn];
      deleteLink(iColumn+numberRows);
      addLink ( iColumn + numberRows, number );
      }
      lastColumn[maximumColumnsExtra_] = numberColumns - 1;
      nextColumn[maximumColumnsExtra_] = 0;
      lastColumn[0] = maximumColumnsExtra_;
      if (numberColumns)
      nextColumn[numberColumns - 1] = maximumColumnsExtra_;
      startColumn[maximumColumnsExtra_] = numberElements;
    }
  }                     /* endswitch */
}

//Does most of factorization
int
00952 CoinFactorization::factor (  )
{
  int * lastColumn = lastColumn_.array();
  int * lastRow = lastRow_.array();
  //sparse
  status_ = factorSparse (  );
  switch ( status_ ) {
  case 0:               //finished
    totalElements_ = 0;
    {
      int * pivotColumn = pivotColumn_.array();
      if ( numberGoodU_ < numberRows_ ) {
      int i, k;
      // Clean out unset nextRow
      int * nextRow = nextRow_.array();
      //int nSing =0;
      k=nextRow[maximumRowsExtra_];
      while (k!=maximumRowsExtra_) {
        int iRow = k;
        k=nextRow[k];
        //nSing++;
        nextRow[iRow]=-1;
      }
      //assert (nSing);
      //printf("%d singularities - good %d rows %d\n",nSing,
      //     numberGoodU_,numberRows_);
      // Now nextRow has -1 or sequence into numberGoodU_;
      int * permuteA = permute_.array();
      for ( i = 0; i < numberRows_; i++ ) {
        int iGood= nextRow[i];
        if (iGood>=0)
          permuteA[iGood]=i;
      }

      // swap arrays
      permute_.swap(nextRow_);
      int * permute = permute_.array();
      for ( i = 0; i < numberRows_; i++ ) {
        lastRow[i] = -1;
      }
      for ( i = 0; i < numberColumns_; i++ ) {
        lastColumn[i] = -1;
      }
      for ( i = 0; i < numberGoodU_; i++ ) {
        int goodRow = permuteA[i];  //valid pivot row
        int goodColumn = pivotColumn[i];

        lastRow[goodRow] = goodColumn;    //will now have -1 or column sequence
        lastColumn[goodColumn] = goodRow; //will now have -1 or row sequence
      }
      nextRow_.conditionalDelete();
      k = 0;
      //copy back and count
      for ( i = 0; i < numberRows_; i++ ) {
        permute[i] = lastRow[i];
        if ( permute[i] < 0 ) {
          //std::cout << i << " " <<permute[i] << std::endl;
        } else {
          k++;
        }
      }
      for ( i = 0; i < numberColumns_; i++ ) {
        pivotColumn[i] = lastColumn[i];
      }
      if ((messageLevel_&4)!=0) 
        std::cout <<"Factorization has "<<numberRows_-k
                <<" singularities"<<std::endl;
      status_ = -1;
      }
    }
    break;
    // dense
  case 2:
    status_=factorDense();
    if(!status_) 
      break;
  default:
    //singular ? or some error
    if ((messageLevel_&4)!=0) 
      std::cout << "Error " << status_ << std::endl;
    break;
  }                     /* endswitch */
  //clean up
  if ( !status_ ) {
    if ( (messageLevel_ & 16)&&numberCompressions_)
      std::cout<<"        Factorization did "<<numberCompressions_
             <<" compressions"<<std::endl;
    if ( numberCompressions_ > 10 ) {
      areaFactor_ *= 1.1;
    }
    numberCompressions_=0;
    cleanup (  );
  }
  return status_;
}


//  pivotRowSingleton.  Does one pivot on Row Singleton in factorization
bool
01051 CoinFactorization::pivotRowSingleton ( int pivotRow,
                          int pivotColumn )
{
  //store pivot columns (so can easily compress)
  CoinBigIndex * startColumnU = startColumnU_.array();
  CoinBigIndex startColumn = startColumnU[pivotColumn];
  int * numberInRow = numberInRow_.array();
  int * numberInColumn = numberInColumn_.array();
  int numberDoColumn = numberInColumn[pivotColumn] - 1;
  CoinBigIndex endColumn = startColumn + numberDoColumn + 1;
  CoinBigIndex pivotRowPosition = startColumn;
  int * indexRowU = indexRowU_.array();
  int iRow = indexRowU[pivotRowPosition];
  CoinBigIndex * startRowU = startRowU_.array();
  int * nextRow = nextRow_.array();
  int * lastRow = lastRow_.array();

  while ( iRow != pivotRow ) {
    pivotRowPosition++;
    iRow = indexRowU[pivotRowPosition];
  }                     /* endwhile */
  assert ( pivotRowPosition < endColumn );
  //store column in L, compress in U and take column out
  CoinBigIndex l = lengthL_;

  if ( l + numberDoColumn > lengthAreaL_ ) {
    //need more memory
    if ((messageLevel_&4)!=0) 
      std::cout << "more memory needed in middle of invert" << std::endl;
    return false;
  }
  CoinBigIndex * startColumnL = startColumnL_.array();
  CoinFactorizationDouble * elementL = elementL_.array();
  int * indexRowL = indexRowL_.array();
  startColumnL[numberGoodL_] = l;   //for luck and first time
  numberGoodL_++;
  startColumnL[numberGoodL_] = l + numberDoColumn;
  lengthL_ += numberDoColumn;
  CoinFactorizationDouble *elementU = elementU_.array();
  CoinFactorizationDouble pivotElement = elementU[pivotRowPosition];
  CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;

  pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
  CoinBigIndex i;

  int * indexColumnU = indexColumnU_.array();
  for ( i = startColumn; i < pivotRowPosition; i++ ) {
    int iRow = indexRowU[i];

    indexRowL[l] = iRow;
    elementL[l] = elementU[i] * pivotMultiplier;
    l++;
    //take out of row list
    CoinBigIndex start = startRowU[iRow];
    int iNumberInRow = numberInRow[iRow];
    CoinBigIndex end = start + iNumberInRow;
    CoinBigIndex where = start;

    while ( indexColumnU[where] != pivotColumn ) {
      where++;
    }                   /* endwhile */
    assert ( where < end );
    indexColumnU[where] = indexColumnU[end - 1];
    iNumberInRow--;
    numberInRow[iRow] = iNumberInRow;
    deleteLink ( iRow );
    addLink ( iRow, iNumberInRow );
  }
  for ( i = pivotRowPosition + 1; i < endColumn; i++ ) {
    int iRow = indexRowU[i];

    indexRowL[l] = iRow;
    elementL[l] = elementU[i] * pivotMultiplier;
    l++;
    //take out of row list
    CoinBigIndex start = startRowU[iRow];
    int iNumberInRow = numberInRow[iRow];
    CoinBigIndex end = start + iNumberInRow;
    CoinBigIndex where = start;

    while ( indexColumnU[where] != pivotColumn ) {
      where++;
    }                   /* endwhile */
    assert ( where < end );
    indexColumnU[where] = indexColumnU[end - 1];
    iNumberInRow--;
    numberInRow[iRow] = iNumberInRow;
    deleteLink ( iRow );
    addLink ( iRow, iNumberInRow );
  }
  numberInColumn[pivotColumn] = 0;
  //modify linked list for pivots
  numberInRow[pivotRow] = 0;
  deleteLink ( pivotRow );
  deleteLink ( pivotColumn + numberRows_ );
  //take out this bit of indexColumnU
  int next = nextRow[pivotRow];
  int last = lastRow[pivotRow];

  nextRow[last] = next;
  lastRow[next] = last;
  lastRow[pivotRow] =-2;
  nextRow[pivotRow] = numberGoodU_; //use for permute
  return true;
}

//  pivotColumnSingleton.  Does one pivot on Column Singleton in factorization
bool
01159 CoinFactorization::pivotColumnSingleton ( int pivotRow,
                             int pivotColumn )
{
  int * numberInRow = numberInRow_.array();
  int * numberInColumn = numberInColumn_.array();
  int * numberInColumnPlus = numberInColumnPlus_.array();
  //store pivot columns (so can easily compress)
  int numberDoRow = numberInRow[pivotRow] - 1;
  CoinBigIndex * startColumnU = startColumnU_.array();
  CoinBigIndex startColumn = startColumnU[pivotColumn];
  int put = 0;
  CoinBigIndex * startRowU = startRowU_.array();
  CoinBigIndex startRow = startRowU[pivotRow];
  CoinBigIndex endRow = startRow + numberDoRow + 1;
  int * indexColumnU = indexColumnU_.array();
  int * indexRowU = indexRowU_.array();
  int * saveColumn = saveColumn_.array();
  CoinBigIndex i;

  for ( i = startRow; i < endRow; i++ ) {
    int iColumn = indexColumnU[i];

    if ( iColumn != pivotColumn ) {
      saveColumn[put++] = iColumn;
    }
  }
  int * nextRow = nextRow_.array();
  int * lastRow = lastRow_.array();
  //take out this bit of indexColumnU
  int next = nextRow[pivotRow];
  int last = lastRow[pivotRow];

  nextRow[last] = next;
  lastRow[next] = last;
  nextRow[pivotRow] = numberGoodU_; //use for permute
  lastRow[pivotRow] =-2; //mark
  //clean up counts
  CoinFactorizationDouble *elementU = elementU_.array();
  CoinFactorizationDouble pivotElement = elementU[startColumn];

  pivotRegion_.array()[numberGoodU_] = 1.0 / pivotElement;
  numberInColumn[pivotColumn] = 0;
  //totalElements_ --;
  //numberInColumnPlus[pivotColumn]++;
  //move pivot row in other columns to safe zone
  for ( i = 0; i < numberDoRow; i++ ) {
    int iColumn = saveColumn[i];

    if ( numberInColumn[iColumn] ) {
      int number = numberInColumn[iColumn] - 1;

      //modify linked list
      deleteLink ( iColumn + numberRows_ );
      addLink ( iColumn + numberRows_, number );
      //move pivot row element
      if ( number ) {
      CoinBigIndex start = startColumnU[iColumn];
      CoinBigIndex pivot = start;
      int iRow = indexRowU[pivot];
      while ( iRow != pivotRow ) {
        pivot++;
        iRow = indexRowU[pivot];
      }
        assert (pivot < startColumnU[iColumn] +
                numberInColumn[iColumn]);
      if ( pivot != start ) {
        //move largest one up
        CoinFactorizationDouble value = elementU[start];

        iRow = indexRowU[start];
        elementU[start] = elementU[pivot];
        indexRowU[start] = indexRowU[pivot];
        elementU[pivot] = elementU[start + 1];
        indexRowU[pivot] = indexRowU[start + 1];
        elementU[start + 1] = value;
        indexRowU[start + 1] = iRow;
      } else {
        //find new largest element
        int iRowSave = indexRowU[start + 1];
        CoinFactorizationDouble valueSave = elementU[start + 1];
        double valueLargest = fabs ( valueSave );
        CoinBigIndex end = start + numberInColumn[iColumn];
        CoinBigIndex largest = start + 1;

        CoinBigIndex k;
        for ( k = start + 2; k < end; k++ ) {
          CoinFactorizationDouble value = elementU[k];
          double valueAbs = fabs ( value );

          if ( valueAbs > valueLargest ) {
            valueLargest = valueAbs;
            largest = k;
          }
        }
        indexRowU[start + 1] = indexRowU[largest];
        elementU[start + 1] = elementU[largest];
        indexRowU[largest] = iRowSave;
        elementU[largest] = valueSave;
      }
      }
      //clean up counts
      numberInColumn[iColumn]--;
      numberInColumnPlus[iColumn]++;
      startColumnU[iColumn]++;
      //totalElements_--;
    }
  }
  //modify linked list for pivots
  deleteLink ( pivotRow );
  deleteLink ( pivotColumn + numberRows_ );
  numberInRow[pivotRow] = 0;
  //put in dummy pivot in L
  CoinBigIndex l = lengthL_;

  CoinBigIndex * startColumnL = startColumnL_.array();
  startColumnL[numberGoodL_] = l;   //for luck and first time
  numberGoodL_++;
  startColumnL[numberGoodL_] = l;
  return true;
}


//  getColumnSpace.  Gets space for one Column with given length
//may have to do compression  (returns true)
//also moves existing vector
bool
01285 CoinFactorization::getColumnSpace ( int iColumn,
                         int extraNeeded )
{
  int * numberInColumn = numberInColumn_.array();
  int * numberInColumnPlus = numberInColumnPlus_.array();
  int * nextColumn = nextColumn_.array();
  int * lastColumn = lastColumn_.array();
  int number = numberInColumnPlus[iColumn] +
    numberInColumn[iColumn];
  CoinBigIndex * startColumnU = startColumnU_.array();
  CoinBigIndex space = lengthAreaU_ - startColumnU[maximumColumnsExtra_];
  CoinFactorizationDouble *elementU = elementU_.array();
  int * indexRowU = indexRowU_.array();

  if ( space < extraNeeded + number + 2 ) {
    //compression
    int iColumn = nextColumn[maximumColumnsExtra_];
    CoinBigIndex put = 0;

    while ( iColumn != maximumColumnsExtra_ ) {
      //move
      CoinBigIndex get;
      CoinBigIndex getEnd;

      if ( startColumnU[iColumn] >= 0 ) {
      get = startColumnU[iColumn]
        - numberInColumnPlus[iColumn];
      getEnd = startColumnU[iColumn] + numberInColumn[iColumn];
      startColumnU[iColumn] = put + numberInColumnPlus[iColumn];
      } else {
      get = -startColumnU[iColumn];
      getEnd = get + numberInColumn[iColumn];
      startColumnU[iColumn] = -put;
      }
      CoinBigIndex i;
      for ( i = get; i < getEnd; i++ ) {
      indexRowU[put] = indexRowU[i];
      elementU[put] = elementU[i];
      put++;
      }
      iColumn = nextColumn[iColumn];
    }                   /* endwhile */
    numberCompressions_++;
    startColumnU[maximumColumnsExtra_] = put;
    space = lengthAreaU_ - put;
    if ( extraNeeded == COIN_INT_MAX >> 1 ) {
      return true;
    }
    if ( space < extraNeeded + number + 2 ) {
      //need more space
      //if we can allocate bigger then do so and copy
      //if not then return so code can start again
      status_ = -99;
      return false;
    }
  }
  CoinBigIndex put = startColumnU[maximumColumnsExtra_];
  int next = nextColumn[iColumn];
  int last = lastColumn[iColumn];

  if ( extraNeeded || next != maximumColumnsExtra_ ) {
    //out
    nextColumn[last] = next;
    lastColumn[next] = last;
    //in at end
    last = lastColumn[maximumColumnsExtra_];
    nextColumn[last] = iColumn;
    lastColumn[maximumColumnsExtra_] = iColumn;
    lastColumn[iColumn] = last;
    nextColumn[iColumn] = maximumColumnsExtra_;
    //move
    CoinBigIndex get = startColumnU[iColumn]
      - numberInColumnPlus[iColumn];

    startColumnU[iColumn] = put + numberInColumnPlus[iColumn];
    if ( number < 50 ) {
      int *indexRow = indexRowU;
      CoinFactorizationDouble *element = elementU;
      int i = 0;

      if ( ( number & 1 ) != 0 ) {
      element[put] = element[get];
      indexRow[put] = indexRow[get];
      i = 1;
      }
      for ( ; i < number; i += 2 ) {
      CoinFactorizationDouble value0 = element[get + i];
      CoinFactorizationDouble value1 = element[get + i + 1];
      int index0 = indexRow[get + i];
      int index1 = indexRow[get + i + 1];

      element[put + i] = value0;
      element[put + i + 1] = value1;
      indexRow[put + i] = index0;
      indexRow[put + i + 1] = index1;
      }
    } else {
      CoinMemcpyN ( &indexRowU[get], number, &indexRowU[put] );
      CoinMemcpyN ( &elementU[get], number, &elementU[put] );
    }
    put += number;
    get += number;
    //add 4 for luck
    startColumnU[maximumColumnsExtra_] = put + extraNeeded + 4;
  } else {
    //take off space
    startColumnU[maximumColumnsExtra_] = startColumnU[last] +
      numberInColumn[last];
  }
  return true;
}

//  getRowSpace.  Gets space for one Row with given length
//may have to do compression  (returns true)
//also moves existing vector
bool
01401 CoinFactorization::getRowSpace ( int iRow,
                      int extraNeeded )
{
  int * numberInRow = numberInRow_.array();
  int number = numberInRow[iRow];
  CoinBigIndex * startRowU = startRowU_.array();
  CoinBigIndex space = lengthAreaU_ - startRowU[maximumRowsExtra_];
  int * nextRow = nextRow_.array();
  int * lastRow = lastRow_.array();
  int * indexColumnU = indexColumnU_.array();

  if ( space < extraNeeded + number + 2 ) {
    //compression
    int iRow = nextRow[maximumRowsExtra_];
    CoinBigIndex put = 0;

    while ( iRow != maximumRowsExtra_ ) {
      //move
      CoinBigIndex get = startRowU[iRow];
      CoinBigIndex getEnd = startRowU[iRow] + numberInRow[iRow];

      startRowU[iRow] = put;
      CoinBigIndex i;
      for ( i = get; i < getEnd; i++ ) {
      indexColumnU[put] = indexColumnU[i];
      put++;
      }
      iRow = nextRow[iRow];
    }                   /* endwhile */
    numberCompressions_++;
    startRowU[maximumRowsExtra_] = put;
    space = lengthAreaU_ - put;
    if ( space < extraNeeded + number + 2 ) {
      //need more space
      //if we can allocate bigger then do so and copy
      //if not then return so code can start again
      status_ = -99;
      return false;;
    }
  }
  CoinBigIndex put = startRowU[maximumRowsExtra_];
  int next = nextRow[iRow];
  int last = lastRow[iRow];

  //out
  nextRow[last] = next;
  lastRow[next] = last;
  //in at end
  last = lastRow[maximumRowsExtra_];
  nextRow[last] = iRow;
  lastRow[maximumRowsExtra_] = iRow;
  lastRow[iRow] = last;
  nextRow[iRow] = maximumRowsExtra_;
  //move
  CoinBigIndex get = startRowU[iRow];

  startRowU[iRow] = put;
  while ( number ) {
    number--;
    indexColumnU[put] = indexColumnU[get];
    put++;
    get++;
  }                     /* endwhile */
  //add 4 for luck
  startRowU[maximumRowsExtra_] = put + extraNeeded + 4;
  return true;
}

#if COIN_ONE_ETA_COPY
/* Reorders U so contiguous and in order (if there is space)
   Returns true if it could */
bool 
CoinFactorization::reorderU()
{
#if 1
  return false;
#else
  if (numberRows_!=numberColumns_)
    return false;
  CoinBigIndex * startColumnU = startColumnU_.array();
  int * numberInColumn = numberInColumn_.array();
  int * numberInColumnPlus = numberInColumnPlus_.array();
  int iColumn;
  CoinBigIndex put = 0;
  for (iColumn =0;iColumn<numberRows_;iColumn++) 
    put += numberInColumnPlus[iColumn];
  CoinBigIndex space = lengthAreaU_ - startColumnU[maximumColumnsExtra_];
  if (space<put) {
    //printf("Space %d out of %d - needed %d\n",
    //   space,lengthAreaU_,put);
    return false;
  }
  int *indexRowU = indexRowU_.array();
  CoinFactorizationDouble *elementU = elementU_.array();
  int * pivotColumn = pivotColumn_.array();
  put = startColumnU[maximumColumnsExtra_];
  for (int jColumn =0;jColumn<numberRows_;jColumn++) {
    iColumn = pivotColumn[jColumn];
    int n = numberInColumnPlus[iColumn];
    CoinBigIndex getEnd = startColumnU[iColumn];
    CoinBigIndex get = getEnd - n;
    startColumnU[iColumn] = put;
    numberInColumn[jColumn]=n;
    CoinBigIndex i;
    for ( i = get; i < getEnd; i++ ) {
      indexRowU[put] = indexRowU[i];
      elementU[put] = elementU[i];
      put++;
    }
  }
  // and pack down
  put = 0;
  for (int jColumn =0;jColumn<numberRows_;jColumn++) {
    iColumn = pivotColumn[jColumn];
    int n = numberInColumn[jColumn];
    CoinBigIndex get = startColumnU[iColumn];
    CoinBigIndex getEnd = get+n;
    CoinBigIndex i;
    for ( i = get; i < getEnd; i++ ) {
      indexRowU[put] = indexRowU[i];
      elementU[put] = elementU[i];
      put++;
    }
  }
  put=0;
  for (iColumn =0;iColumn<numberRows_;iColumn++) { 
    int n = numberInColumn[iColumn];
    startColumnU[iColumn]=put;
    put += n;
    //numberInColumnPlus[iColumn]=n;
    //numberInColumn[iColumn]=0; // necessary?
    //pivotColumn[iColumn]=iColumn;
  }
#if 0
  // reset nextColumn - probably not necessary
  int * nextColumn = nextColumn_.array();
  nextColumn[maximumColumnsExtra_]=0;
  nextColumn[numberRows_-1] = maximumColumnsExtra_;
  for (iColumn=0;iColumn<numberRows_-1;iColumn++)
    nextColumn[iColumn]=iColumn+1;
  // swap arrays
  numberInColumn_.swap(numberInColumnPlus_);
#endif
  //return false;
  return true;
#endif
}
#endif
//  cleanup.  End of factorization
void
01551 CoinFactorization::cleanup (  )
{
#if COIN_ONE_ETA_COPY
  bool compressDone = reorderU();
  if (!compressDone) {
    getColumnSpace ( 0, COIN_INT_MAX >> 1 );    //compress
    // swap arrays
    numberInColumn_.swap(numberInColumnPlus_);
  }
#else
  getColumnSpace ( 0, COIN_INT_MAX >> 1 );      //compress
  // swap arrays
    numberInColumn_.swap(numberInColumnPlus_);
#endif
  CoinBigIndex * startColumnU = startColumnU_.array();
  CoinBigIndex lastU = startColumnU[maximumColumnsExtra_];

  //free some memory here
  saveColumn_.conditionalDelete();
  markRow_.conditionalDelete() ;
  //firstCount_.conditionalDelete() ;
  nextCount_.conditionalDelete() ;
  lastCount_.conditionalDelete() ;
  int * numberInRow = numberInRow_.array();
  int * numberInColumn = numberInColumn_.array();
  int * numberInColumnPlus = numberInColumnPlus_.array();
  //make column starts OK
  //for best cache behavior get in order (last pivot at bottom of space)
  //that will need thinking about
  //use nextRow for permutation  (as that is what it is)
  int i;

#ifndef NDEBUG
  {
    if (numberGoodU_<numberRows_)
      abort();
    char * mark = new char[numberRows_];
    memset(mark,0,numberRows_);
    int * array;
    array = nextRow_.array();
    for ( i = 0; i < numberRows_; i++ ) {
      int k = array[i];
      if(k<0||k>=numberRows_)
      printf("Bad a %d %d\n",i,k);
      assert(k>=0&&k<numberRows_);
      if(mark[k]==1)
      printf("Bad a %d %d\n",i,k);
      mark[k]=1;
    }
    for ( i = 0; i < numberRows_; i++ ) {
      assert(mark[i]==1);
      if(mark[i]!=1)
      printf("Bad b %d\n",i);
    }
    delete [] mark;
  }
#endif
  // swap arrays
  permute_.swap(nextRow_);
  //safety feature
  int * permute = permute_.array();
  permute[numberRows_] = 0;
  permuteBack_.conditionalNew( maximumRowsExtra_ + 1);
  int * permuteBack = permuteBack_.array();
#ifdef ZEROFAULT
  memset(permuteBack_.array(),'w',(maximumRowsExtra_+1)*sizeof(int));
#endif
  for ( i = 0; i < numberRows_; i++ ) {
    int iRow = permute[i];

    permuteBack[iRow] = i;
  }
  //redo nextRow_

#ifndef NDEBUG
  for ( i = 0; i < numberRows_; i++ ) {
    assert (permute[i]>=0&&permute[i]<numberRows_);
    assert (permuteBack[i]>=0&&permuteBack[i]<numberRows_);
  }
#endif
#if COIN_ONE_ETA_COPY
  if (!compressDone) {
#endif
    // Redo total elements
    totalElements_=0;
    for ( i = 0; i < numberColumns_; i++ ) {
      int number = numberInColumn[i];
      totalElements_ += number;
      startColumnU[i] -= number;
    }
#if COIN_ONE_ETA_COPY
  }
#endif
  int numberU = 0;

  pivotColumnBack_.conditionalNew( maximumRowsExtra_ + 1);
#ifdef ZEROFAULT
  memset(pivotColumnBack(),'q',(maximumRowsExtra_+1)*sizeof(int));
#endif
  const int * pivotColumn = pivotColumn_.array();
  int * pivotColumnB = pivotColumnBack();
  int *indexColumnU = indexColumnU_.array();
  CoinBigIndex *startColumn = startColumnU;
  int *indexRowU = indexRowU_.array();
  CoinFactorizationDouble *elementU = elementU_.array();
#if COIN_ONE_ETA_COPY
  if (!compressDone) {
#endif
    for ( i = 0; i < numberColumns_; i++ ) {
      int iColumn = pivotColumn[i];
      
      pivotColumnB[iColumn] = i;
      if ( iColumn >= 0 ) {
      //wanted
      if ( numberU != iColumn ) {
        numberInColumnPlus[iColumn] = numberU;
      } else {
        numberInColumnPlus[iColumn] = -1; //already in correct place
      }
      numberU++;
      }
    }
    for ( i = 0; i < numberColumns_; i++ ) {
      int number = numberInColumn[i];     //always 0?
      int where = numberInColumnPlus[i];
      
      numberInColumnPlus[i] = -1;
      CoinBigIndex start = startColumnU[i];
      
      while ( where >= 0 ) {
      //put where it should be
      int numberNext = numberInColumn[where];   //always 0?
      int whereNext = numberInColumnPlus[where];
      CoinBigIndex startNext = startColumnU[where];
      
      numberInColumn[where] = number;
      numberInColumnPlus[where] = -1;
      startColumnU[where] = start;
      number = numberNext;
      where = whereNext;
      start = startNext;
      }                       /* endwhile */
    }
    //sort - using indexColumn
    CoinFillN ( indexColumnU_.array(), lastU, -1 );
    CoinBigIndex k = 0;
    
    for ( i = numberSlacks_; i < numberRows_; i++ ) {
      CoinBigIndex start = startColumn[i];
      CoinBigIndex end = start + numberInColumn[i];
      
      CoinBigIndex j;
      for ( j = start; j < end; j++ ) {
      indexColumnU[j] = k++;
      }
    }
    for ( i = numberSlacks_; i < numberRows_; i++ ) {
      CoinBigIndex start = startColumn[i];
      CoinBigIndex end = start + numberInColumn[i];
      
      CoinBigIndex j;
      for ( j = start; j < end; j++ ) {
      CoinBigIndex k = indexColumnU[j];
      int iRow = indexRowU[j];
      CoinFactorizationDouble element = elementU[j];
      
      while ( k != -1 ) {
        CoinBigIndex kNext = indexColumnU[k];
        int iRowNext = indexRowU[k];
        CoinFactorizationDouble elementNext = elementU[k];
        
        indexColumnU[k] = -1;
        indexRowU[k] = iRow;
        elementU[k] = element;
        k = kNext;
        iRow = iRowNext;
        element = elementNext;
      }                       /* endwhile */
      }
    }
    CoinZeroN ( startColumnU, numberSlacks_ );
    k = 0;
    for ( i = numberSlacks_; i < numberRows_; i++ ) {
      startColumnU[i] = k;
      k += numberInColumn[i];
    }
    maximumU_=k;
#if COIN_ONE_ETA_COPY
  } else {
    // U already OK
    for ( i = 0; i < numberColumns_; i++ ) {
      int iColumn = pivotColumn[i];
      pivotColumnB[iColumn] = i;
    }
    maximumU_=startColumnU[numberRows_-1]+
      numberInColumn[numberRows_-1];
    numberU=numberRows_;
  }
#endif
  if ( (messageLevel_ & 8)) {
    std::cout<<"        length of U "<<totalElements_<<", length of L "<<lengthL_;
    if (numberDense_)
      std::cout<<" plus "<<numberDense_*numberDense_<<" from "<<numberDense_<<" dense rows";
    std::cout<<std::endl;
  }
  // and add L and dense
  totalElements_ += numberDense_*numberDense_+lengthL_;
  int * nextColumn = nextColumn_.array();
  int * lastColumn = lastColumn_.array();
  // See whether to have extra copy of R
  if (maximumU_>10*numberRows_||numberRows_<200) {
    // NO
    numberInColumnPlus_.conditionalDelete() ;
  } else {
    for ( i = 0; i < numberColumns_; i++ ) {
      lastColumn[i] = i - 1;
      nextColumn[i] = i + 1;
      numberInColumnPlus[i]=0;
    }
    nextColumn[numberColumns_ - 1] = maximumColumnsExtra_;
    lastColumn[maximumColumnsExtra_] = numberColumns_ - 1;
    nextColumn[maximumColumnsExtra_] = 0;
    lastColumn[0] = maximumColumnsExtra_;
  }
  numberU_ = numberU;
  numberGoodU_ = numberU;
  numberL_ = numberGoodL_;
#if COIN_DEBUG
  for ( i = 0; i < numberRows_; i++ ) {
    if ( permute[i] < 0 ) {
      std::cout << i << std::endl;
      abort (  );
    }
  }
#endif
  CoinFactorizationDouble * pivotRegion = pivotRegion_.array();
  for ( i = numberSlacks_; i < numberU; i++ ) {
    CoinBigIndex start = startColumnU[i];
    CoinBigIndex end = start + numberInColumn[i];

    totalElements_ += numberInColumn[i];
    CoinBigIndex j;

    for ( j = start; j < end; j++ ) {
      int iRow = indexRowU[j];
      iRow = permute[iRow];
      indexRowU[j] = iRow;
      numberInRow[iRow]++;
    }
  }
#if COIN_ONE_ETA_COPY
  if (numberRows_>=COIN_ONE_ETA_COPY) {
#endif
    //space for cross reference
    convertRowToColumnU_.conditionalNew( lengthAreaU_ );
    CoinBigIndex *convertRowToColumn = convertRowToColumnU_.array();
    CoinBigIndex j = 0;
    CoinBigIndex *startRow = startRowU_.array();
    
    int iRow;
    for ( iRow = 0; iRow < numberRows_; iRow++ ) {
      startRow[iRow] = j;
      j += numberInRow[iRow];
    }
    CoinBigIndex numberInU = j;
    
    CoinZeroN ( numberInRow_.array(), numberRows_ );
    
    for ( i = numberSlacks_; i < numberRows_; i++ ) {
      CoinBigIndex start = startColumnU[i];
      CoinBigIndex end = start + numberInColumn[i];
      
      CoinFactorizationDouble pivotValue = pivotRegion[i];

      CoinBigIndex j;
      for ( j = start; j < end; j++ ) {
      int iRow = indexRowU[j];
      int iLook = numberInRow[iRow];
      
      numberInRow[iRow] = iLook + 1;
      CoinBigIndex k = startRow[iRow] + iLook;
      
      indexColumnU[k] = i;
      convertRowToColumn[k] = j;
      //multiply by pivot
      elementU[j] *= pivotValue;
      }
    }
    int * nextRow = nextRow_.array();
    int * lastRow = lastRow_.array();
    for ( j = 0; j < numberRows_; j++ ) {
      lastRow[j] = j - 1;
      nextRow[j] = j + 1;
    }
    nextRow[numberRows_ - 1] = maximumRowsExtra_;
    lastRow[maximumRowsExtra_] = numberRows_ - 1;
    nextRow[maximumRowsExtra_] = 0;
    lastRow[0] = maximumRowsExtra_;
    startRow[maximumRowsExtra_] = numberInU;
#if COIN_ONE_ETA_COPY
  } else {
    // no row copy
    for ( i = numberSlacks_; i < numberU; i++ ) {
      CoinBigIndex start = startColumnU[i];
      CoinBigIndex end = start + numberInColumn[i];
      
      CoinBigIndex j;
      CoinFactorizationDouble pivotValue = pivotRegion[i];
      
      for ( j = start; j < end; j++ ) {
      //multiply by pivot
      elementU[j] *= pivotValue;
      }
    }
  }
#endif

  int firstReal = numberRows_;

  CoinBigIndex * startColumnL = startColumnL_.array();
  int * indexRowL = indexRowL_.array();
  for ( i = numberRows_ - 1; i >= 0; i-- ) {
    CoinBigIndex start = startColumnL[i];
    CoinBigIndex end = startColumnL[i + 1];

    totalElements_ += end - start;
    if ( end > start ) {
      firstReal = i;
      CoinBigIndex j;
      for ( j = start; j < end; j++ ) {
      int iRow = indexRowL[j];
      iRow = permute[iRow];
      assert (iRow>firstReal);
      indexRowL[j] = iRow;
      }
    }
  }
  baseL_ = firstReal;
  numberL_ -= firstReal;
  factorElements_ = totalElements_;
  //can delete pivotRowL_ as not used
  pivotRowL_.conditionalDelete() ;
  //use L for R if room
  CoinBigIndex space = lengthAreaL_ - lengthL_;
  CoinBigIndex spaceUsed = lengthL_ + lengthU_;

  int needed = ( spaceUsed + numberRows_ - 1 ) / numberRows_;

  needed = needed * 2 * maximumPivots_;
  if ( needed < 2 * numberRows_ ) {
    needed = 2 * numberRows_;
  }
  if (numberInColumnPlus_.array()) {
    // Need double the space for R
    space = space/2;
    startColumnR_.conditionalNew( maximumPivots_ + 1 + maximumColumnsExtra_ + 1 );
    CoinBigIndex * startR = startColumnR_.array() + maximumPivots_+1;
    CoinZeroN (startR,(maximumColumnsExtra_+1));
  } else {
    startColumnR_.conditionalNew(maximumPivots_ + 1 );
  }
#ifdef ZEROFAULT
  memset(startColumnR_.array(),'z',(maximumPivots_ + 1)*sizeof(CoinBigIndex));
#endif
  if ( space >= needed ) {
    lengthR_ = 0;
    lengthAreaR_ = space;
    elementR_ = elementL_.array() + lengthL_;
    indexRowR_ = indexRowL_.array() + lengthL_;
  } else {
    lengthR_ = 0;
    lengthAreaR_ = space;
    elementR_ = elementL_.array() + lengthL_;
    indexRowR_ = indexRowL_.array() + lengthL_;
    if ((messageLevel_&4))
      std::cout<<"Factorization may need some increasing area space"
             <<std::endl;
    if ( areaFactor_ ) {
      areaFactor_ *= 1.1;
    } else {
      areaFactor_ = 1.1;
    }
  }
  numberR_ = 0;
}
// Returns areaFactor but adjusted for dense
double 
01938 CoinFactorization::adjustedAreaFactor() const
{
  double factor = areaFactor_;
  if (numberDense_&&areaFactor_>1.0) {
    double dense = numberDense_;
    dense *= dense;
    double withoutDense = totalElements_ - dense +1.0;
    factor *= 1.0 +dense/withoutDense;
  }
  return factor;
}

//  checkConsistency.  Checks that row and column copies look OK
void
01952 CoinFactorization::checkConsistency (  )
{
  bool bad = false;

  int iRow;
  CoinBigIndex * startRowU = startRowU_.array();
  int * numberInRow = numberInRow_.array();
  int * numberInColumn = numberInColumn_.array();
  int * indexColumnU = indexColumnU_.array();
  int * indexRowU = indexRowU_.array();
  CoinBigIndex * startColumnU = startColumnU_.array();
  for ( iRow = 0; iRow < numberRows_; iRow++ ) {
    if ( numberInRow[iRow] ) {
      CoinBigIndex startRow = startRowU[iRow];
      CoinBigIndex endRow = startRow + numberInRow[iRow];

      CoinBigIndex j;
      for ( j = startRow; j < endRow; j++ ) {
      int iColumn = indexColumnU[j];
      CoinBigIndex startColumn = startColumnU[iColumn];
      CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
      bool found = false;

      CoinBigIndex k;
      for ( k = startColumn; k < endColumn; k++ ) {
        if ( indexRowU[k] == iRow ) {
          found = true;
          break;
        }
      }
      if ( !found ) {
        bad = true;
        std::cout << "row " << iRow << " column " << iColumn << " Rows" << std::endl;
      }
      }
    }
  }
  int iColumn;
  for ( iColumn = 0; iColumn < numberColumns_; iColumn++ ) {
    if ( numberInColumn[iColumn] ) {
      CoinBigIndex startColumn = startColumnU[iColumn];
      CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];

      CoinBigIndex j;
      for ( j = startColumn; j < endColumn; j++ ) {
      int iRow = indexRowU[j];
      CoinBigIndex startRow = startRowU[iRow];
      CoinBigIndex endRow = startRow + numberInRow[iRow];
      bool found = false;

      CoinBigIndex k;
      for (  k = startRow; k < endRow; k++ ) {
        if ( indexColumnU[k] == iColumn ) {
          found = true;
          break;
        }
      }
      if ( !found ) {
        bad = true;
        std::cout << "row " << iRow << " column " << iColumn << " Columns" <<
          std::endl;
      }
      }
    }
  }
  if ( bad ) {
    abort (  );
  }
}
  //  pivotOneOtherRow.  When just one other row so faster
bool 
02023 CoinFactorization::pivotOneOtherRow ( int pivotRow,
                                 int pivotColumn )
{
  int * numberInRow = numberInRow_.array();
  int * numberInColumn = numberInColumn_.array();
  int * numberInColumnPlus = numberInColumnPlus_.array();
  int numberInPivotRow = numberInRow[pivotRow] - 1;
  CoinBigIndex * startRowU = startRowU_.array();
  CoinBigIndex * startColumnU = startColumnU_.array();
  CoinBigIndex startColumn = startColumnU[pivotColumn];
  CoinBigIndex startRow = startRowU[pivotRow];
  CoinBigIndex endRow = startRow + numberInPivotRow + 1;

  //take out this bit of indexColumnU
  int * nextRow = nextRow_.array();
  int * lastRow = lastRow_.array();
  int next = nextRow[pivotRow];
  int last = lastRow[pivotRow];

  nextRow[last] = next;
  lastRow[next] = last;
  nextRow[pivotRow] = numberGoodU_; //use for permute
  lastRow[pivotRow] = -2;
  numberInRow[pivotRow] = 0;
  //store column in L, compress in U and take column out
  CoinBigIndex l = lengthL_;

  if ( l + 1 > lengthAreaL_ ) {
    //need more memory
    if ((messageLevel_&4)!=0) 
      std::cout << "more memory needed in middle of invert" << std::endl;
    return false;
  }
  //l+=currentAreaL_->elementByColumn-elementL_;
  //CoinBigIndex lSave=l;
  CoinBigIndex * startColumnL = startColumnL_.array();
  CoinFactorizationDouble * elementL = elementL_.array();
  int * indexRowL = indexRowL_.array();
  startColumnL[numberGoodL_] = l;   //for luck and first time
  numberGoodL_++;
  startColumnL[numberGoodL_] = l + 1;
  lengthL_++;
  CoinFactorizationDouble pivotElement;
  CoinFactorizationDouble otherMultiplier;
  int otherRow;
  int * saveColumn = saveColumn_.array();
  CoinFactorizationDouble *elementU = elementU_.array();
  int * indexRowU = indexRowU_.array();

  if ( indexRowU[startColumn] == pivotRow ) {
    pivotElement = elementU[startColumn];
    otherMultiplier = elementU[startColumn + 1];
    otherRow = indexRowU[startColumn + 1];
  } else {
    pivotElement = elementU[startColumn + 1];
    otherMultiplier = elementU[startColumn];
    otherRow = indexRowU[startColumn];
  }
  int numberSave = numberInRow[otherRow];
  CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;

  CoinFactorizationDouble * pivotRegion = pivotRegion_.array();
  pivotRegion[numberGoodU_] = pivotMultiplier;
  numberInColumn[pivotColumn] = 0;
  otherMultiplier = otherMultiplier * pivotMultiplier;
  indexRowL[l] = otherRow;
  elementL[l] = otherMultiplier;
  //take out of row list
  CoinBigIndex start = startRowU[otherRow];
  CoinBigIndex end = start + numberSave;
  CoinBigIndex where = start;
  int * indexColumnU = indexColumnU_.array();

  while ( indexColumnU[where] != pivotColumn ) {
    where++;
  }                     /* endwhile */
  assert ( where < end );
  end--;
  indexColumnU[where] = indexColumnU[end];
  int numberAdded = 0;
  int numberDeleted = 0;

  //pack down and move to work
  int j;
  const int * nextCount = nextCount_.array();
  int * nextColumn = nextColumn_.array();

  for ( j = startRow; j < endRow; j++ ) {
    int iColumn = indexColumnU[j];

    if ( iColumn != pivotColumn ) {
      CoinBigIndex startColumn = startColumnU[iColumn];
      CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
      int iRow = indexRowU[startColumn];
      CoinFactorizationDouble value = elementU[startColumn];
      double largest;
      bool foundOther = false;

      //leave room for pivot
      CoinBigIndex put = startColumn + 1;
      CoinBigIndex positionLargest = -1;
      CoinFactorizationDouble thisPivotValue = 0.0;
      CoinFactorizationDouble otherElement = 0.0;
      CoinFactorizationDouble nextValue = elementU[put];;
      int nextIRow = indexRowU[put];

      //compress column and find largest not updated
      if ( iRow != pivotRow ) {
      if ( iRow != otherRow ) {
        largest = fabs ( value );
        elementU[put] = value;
        indexRowU[put] = iRow;
        positionLargest = put;
        put++;
        CoinBigIndex i;
        for ( i = startColumn + 1; i < endColumn; i++ ) {
          iRow = nextIRow;
          value = nextValue;
#ifdef ZEROFAULT
          // doesn't matter reading uninitialized but annoys checking
          if ( i + 1 < endColumn ) {
#endif
            nextIRow = indexRowU[i + 1];
            nextValue = elementU[i + 1];
#ifdef ZEROFAULT
          }
#endif
          if ( iRow != pivotRow ) {
            if ( iRow != otherRow ) {
            //keep
            indexRowU[put] = iRow;
            elementU[put] = value;;
            put++;
            } else {
            otherElement = value;
            foundOther = true;
            }
          } else {
            thisPivotValue = value;
          }
        }
      } else {
        otherElement = value;
        foundOther = true;
        //need to find largest
        largest = 0.0;
        CoinBigIndex i;
        for ( i = startColumn + 1; i < endColumn; i++ ) {
          iRow = nextIRow;
          value = nextValue;
#ifdef ZEROFAULT
          // doesn't matter reading uninitialized but annoys checking
          if ( i + 1 < endColumn ) {
#endif
            nextIRow = indexRowU[i + 1];
            nextValue = elementU[i + 1];
#ifdef ZEROFAULT
          }
#endif
          if ( iRow != pivotRow ) {
            //keep
            indexRowU[put] = iRow;
            elementU[put] = value;;
            double absValue = fabs ( value );

            if ( absValue > largest ) {
            largest = absValue;
            positionLargest = put;
            }
            put++;
          } else {
            thisPivotValue = value;
          }
        }
      }
      } else {
      //need to find largest
      largest = 0.0;
      thisPivotValue = value;
      CoinBigIndex i;
      for ( i = startColumn + 1; i < endColumn; i++ ) {
        iRow = nextIRow;
        value = nextValue;
#ifdef ZEROFAULT
        // doesn't matter reading uninitialized but annoys checking
        if ( i + 1 < endColumn ) {
#endif
          nextIRow = indexRowU[i + 1];
          nextValue = elementU[i + 1];
#ifdef ZEROFAULT
        }
#endif
        if ( iRow != otherRow ) {
          //keep
          indexRowU[put] = iRow;
          elementU[put] = value;;
          double absValue = fabs ( value );

          if ( absValue > largest ) {
            largest = absValue;
            positionLargest = put;
          }
          put++;
        } else {
          otherElement = value;
          foundOther = true;
        }
      }
      }
      //slot in pivot
      elementU[startColumn] = thisPivotValue;
      indexRowU[startColumn] = pivotRow;
      //clean up counts
      startColumn++;
      numberInColumn[iColumn] = put - startColumn;
      numberInColumnPlus[iColumn]++;
      startColumnU[iColumn]++;
      otherElement = otherElement - thisPivotValue * otherMultiplier;
      double absValue = fabs ( otherElement );

      if ( absValue > zeroTolerance_ ) {
      if ( !foundOther ) {
        //have we space
        saveColumn[numberAdded++] = iColumn;
        int next = nextColumn[iColumn];
        CoinBigIndex space;

        space = startColumnU[next] - put - numberInColumnPlus[next];
        if ( space <= 0 ) {
          //getColumnSpace also moves fixed part
          int number = numberInColumn[iColumn];

          if ( !getColumnSpace ( iColumn, number + 1 ) ) {
            return false;
          }
          //redo starts
          positionLargest =
            positionLargest + startColumnU[iColumn] - startColumn;
          startColumn = startColumnU[iColumn];
          put = startColumn + number;
        }
      }
      elementU[put] = otherElement;
      indexRowU[put] = otherRow;
      if ( absValue > largest ) {
        largest = absValue;
        positionLargest = put;
      }
      put++;
      } else {
      if ( foundOther ) {
        numberDeleted++;
        //take out of row list
        CoinBigIndex where = start;

        while ( indexColumnU[where] != iColumn ) {
          where++;
        }               /* endwhile */
        assert ( where < end );
        end--;
        indexColumnU[where] = indexColumnU[end];
      }
      }
      numberInColumn[iColumn] = put - startColumn;
      //move largest
      if ( positionLargest >= 0 ) {
      value = elementU[positionLargest];
      iRow = indexRowU[positionLargest];
      elementU[positionLargest] = elementU[startColumn];
      indexRowU[positionLargest] = indexRowU[startColumn];
      elementU[startColumn] = value;
      indexRowU[startColumn] = iRow;
      }
      //linked list for column
      if ( nextCount[iColumn + numberRows_] != -2 ) {
      //modify linked list
      deleteLink ( iColumn + numberRows_ );
      addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
      }
    }
  }
  //get space for row list
  next = nextRow[otherRow];
  CoinBigIndex space;

  space = startRowU[next] - end;
  totalElements_ += numberAdded - numberDeleted;
  int number = numberAdded + ( end - start );

  if ( space < numberAdded ) {
    numberInRow[otherRow] = end - start;
    if ( !getRowSpace ( otherRow, number ) ) {
      return false;
    }
    end = startRowU[otherRow] + end - start;
  }
  // do linked lists and update counts
  numberInRow[otherRow] = number;
  if ( number != numberSave ) {
    deleteLink ( otherRow );
    addLink ( otherRow, number );
  }
  for ( j = 0; j < numberAdded; j++ ) {
    indexColumnU[end++] = saveColumn[j];
  }
  //modify linked list for pivots
  deleteLink ( pivotRow );
  deleteLink ( pivotColumn + numberRows_ );
  return true;
}
void 
CoinFactorization::setPersistenceFlag(int flag)
{ 
  persistenceFlag_=flag;
  workArea_.setPersistence(flag,maximumRowsExtra_+1);
  workArea2_.setPersistence(flag,maximumRowsExtra_+1);
  pivotColumn_.setPersistence(flag,maximumColumnsExtra_+1);
  permute_.setPersistence(flag,maximumRowsExtra_+1);
  pivotColumnBack_.setPersistence(flag,maximumRowsExtra_+1);
  permuteBack_.setPersistence(flag,maximumRowsExtra_+1);
  nextRow_.setPersistence(flag,maximumRowsExtra_+1);
  startRowU_.setPersistence(flag,maximumRowsExtra_+1);
  numberInRow_.setPersistence(flag,maximumRowsExtra_+1);
  numberInColumn_.setPersistence(flag,maximumColumnsExtra_+1);
  numberInColumnPlus_.setPersistence(flag,maximumColumnsExtra_+1);
  firstCount_.setPersistence(flag,CoinMax(biggerDimension_+2,maximumRowsExtra_+1));
  nextCount_.setPersistence(flag,numberRows_+numberColumns_);
  lastCount_.setPersistence(flag,numberRows_+numberColumns_);
  nextColumn_.setPersistence(flag,maximumColumnsExtra_+1);
  lastColumn_.setPersistence(flag,maximumColumnsExtra_+1);
  lastRow_.setPersistence(flag,maximumRowsExtra_+1);
  markRow_.setPersistence(flag,numberRows_);
  saveColumn_.setPersistence(flag,numberColumns_);
  indexColumnU_.setPersistence(flag,lengthAreaU_);
  pivotRowL_.setPersistence(flag,numberRows_+1);
  pivotRegion_.setPersistence(flag,maximumRowsExtra_+1);
  elementU_.setPersistence(flag,lengthAreaU_);
  indexRowU_.setPersistence(flag,lengthAreaU_);
  startColumnU_.setPersistence(flag,maximumColumnsExtra_+1);
  convertRowToColumnU_.setPersistence( flag, lengthAreaU_ );
  elementL_.setPersistence( flag, lengthAreaL_ );
  indexRowL_.setPersistence( flag, lengthAreaL_ );
  startColumnL_.setPersistence( flag, numberRows_ + 1 );
  startColumnR_.setPersistence( flag,  maximumPivots_ + 1 + maximumColumnsExtra_ + 1 );
  startRowL_.setPersistence( flag,0 );
  indexColumnL_.setPersistence( flag, 0 );
  elementByRowL_.setPersistence( flag, 0 );
  sparse_.setPersistence( flag, 0 );
}
// Delete all stuff
void 
02374 CoinFactorization::almostDestructor()
{
  gutsOfDestructor(1);
}

Generated by  Doxygen 1.6.0   Back to index