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

CoinFactorization3.cpp

/* $Id: CoinFactorization3.cpp 1237 2009-12-05 22:43:00Z bjarni $ */
// 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 <cstdio>

#include "CoinFactorization.hpp"
#include "CoinIndexedVector.hpp"
#include "CoinHelperFunctions.hpp"
#include <stdio.h>
#include <iostream>
#if DENSE_CODE==1 
// using simple lapack interface
extern "C" 
{
  /** LAPACK Fortran subroutine DGETRS. */
  void F77_FUNC(dgetrs,DGETRS)(char *trans, cipfint *n,
                               cipfint *nrhs, const double *A, cipfint *ldA,
                               cipfint * ipiv, double *B, cipfint *ldB, ipfint *info,
                         int trans_len);
}
#endif
// For semi-sparse
#define BITS_PER_CHECK 8
#define CHECK_SHIFT 3
typedef unsigned char CoinCheckZero;

//:class CoinFactorization.  Deals with Factorization and Updates

/* Updates one column (FTRAN) from region2 and permutes.
   region1 starts as zero
   Note - if regionSparse2 packed on input - will be packed on output
   - returns un-permuted result in region2 and region1 is zero */
00042 int CoinFactorization::updateColumn ( CoinIndexedVector * regionSparse,
                              CoinIndexedVector * regionSparse2,
                              bool noPermute) 
  const
{
  //permute and move indices into index array
  int * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
  int numberNonZero;
  const int *permute = permute_.array();
  double * COIN_RESTRICT region = regionSparse->denseVector();

#ifndef CLP_FACTORIZATION
  if (!noPermute) {
#endif
    numberNonZero = regionSparse2->getNumElements();
    int * COIN_RESTRICT index = regionSparse2->getIndices();
    double * COIN_RESTRICT array = regionSparse2->denseVector();
#ifndef CLP_FACTORIZATION
    bool packed = regionSparse2->packedMode();
    if (packed) {
      for (int j = 0; j < numberNonZero; j ++ ) {
      int iRow = index[j];
      double value = array[j];
      array[j]=0.0;
      iRow = permute[iRow];
      region[iRow] = value;
      regionIndex[j] = iRow;
      }
    } else {
#else
      assert (!regionSparse2->packedMode());
#endif
      for (int j = 0; j < numberNonZero; j ++ ) {
      int iRow = index[j];
      double value = array[iRow];
      array[iRow]=0.0;
      iRow = permute[iRow];
      region[iRow] = value;
      regionIndex[j] = iRow;
      }
#ifndef CLP_FACTORIZATION
    }
#endif
    regionSparse->setNumElements ( numberNonZero );
#ifndef CLP_FACTORIZATION
  } else {
    numberNonZero = regionSparse->getNumElements();
  }
#endif
  if (collectStatistics_) {
    numberFtranCounts_++;
    ftranCountInput_ += numberNonZero;
  }
    
  //  ******* L
  updateColumnL ( regionSparse, regionIndex );
  if (collectStatistics_) 
    ftranCountAfterL_ += regionSparse->getNumElements();
  //permute extra
  //row bits here
  updateColumnR ( regionSparse );
  if (collectStatistics_) 
    ftranCountAfterR_ += regionSparse->getNumElements();
  
  //update counts
  //  ******* U
  updateColumnU ( regionSparse, regionIndex);
  if (!doForrestTomlin_) {
    // Do PFI after everything else
    updateColumnPFI(regionSparse);
  }
  if (!noPermute) {
    permuteBack(regionSparse,regionSparse2);
    return regionSparse2->getNumElements (  );
  } else {
    return regionSparse->getNumElements (  );
  }
}
// Permutes back at end of updateColumn
void 
00122 CoinFactorization::permuteBack ( CoinIndexedVector * regionSparse, 
                         CoinIndexedVector * outVector) const
{
  // permute back
  int oldNumber = regionSparse->getNumElements();
  const int * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
  double * COIN_RESTRICT region = regionSparse->denseVector();
  int * COIN_RESTRICT outIndex = outVector->getIndices (  );
  double * COIN_RESTRICT out = outVector->denseVector();
  const int * COIN_RESTRICT permuteBack = pivotColumnBack();
  int number=0;
  if (outVector->packedMode()) {
    for (int j = 0; j < oldNumber; j ++ ) {
      int iRow = regionIndex[j];
      double value = region[iRow];
      region[iRow]=0.0;
      if (fabs(value)>zeroTolerance_) {
      iRow = permuteBack[iRow];
      outIndex[number]=iRow;
      out[number++] = value;
      }
    }
  } else {
    for (int j = 0; j < oldNumber; j ++ ) {
      int iRow = regionIndex[j];
      double value = region[iRow];
      region[iRow]=0.0;
      if (fabs(value)>zeroTolerance_) {
      iRow = permuteBack[iRow];
      outIndex[number++]=iRow;
      out[iRow] = value;
      }
    }
  }
  outVector->setNumElements(number);
  regionSparse->setNumElements(0);
}
//  updateColumnL.  Updates part of column (FTRANL)
void
CoinFactorization::updateColumnL ( CoinIndexedVector * regionSparse,
                                 int * COIN_RESTRICT regionIndex) const
{
  if (numberL_) {
    int number = regionSparse->getNumElements (  );
    int goSparse;
    // Guess at number at end
    if (sparseThreshold_>0) {
      if (ftranAverageAfterL_) {
      int newNumber = static_cast<int> (number*ftranAverageAfterL_);
      if (newNumber< sparseThreshold_&&(numberL_<<2)>newNumber)
        goSparse = 2;
      else if (newNumber< sparseThreshold2_&&(numberL_<<1)>newNumber)
        goSparse = 1;
      else
        goSparse = 0;
      } else {
      if (number<sparseThreshold_&&(numberL_<<2)>number) 
        goSparse = 2;
      else
        goSparse = 0;
      }
    } else {
      goSparse=0;
    }
    switch (goSparse) {
    case 0: // densish
      updateColumnLDensish(regionSparse,regionIndex);
      break;
    case 1: // middling
      updateColumnLSparsish(regionSparse,regionIndex);
      break;
    case 2: // sparse
      updateColumnLSparse(regionSparse,regionIndex);
      break;
    }
  }
#ifdef DENSE_CODE
  if (numberDense_) {
    //take off list
    int lastSparse = numberRows_-numberDense_;
    int number = regionSparse->getNumElements();
    double * COIN_RESTRICT region = regionSparse->denseVector (  );
    int i=0;
    bool doDense=false;
    while (i<number) {
      int iRow = regionIndex[i];
      if (iRow>=lastSparse) {
      doDense=true;
      regionIndex[i] = regionIndex[--number];
      } else {
      i++;
      }
    }
    if (doDense) {
      char trans = 'N';
      int ione=1;
      int info;
      F77_FUNC(dgetrs,DGETRS)(&trans,&numberDense_,&ione,denseArea_,&numberDense_,
                        densePermute_,region+lastSparse,&numberDense_,&info,1);
      for (int i=lastSparse;i<numberRows_;i++) {
      double value = region[i];
      if (value) {
        if (fabs(value)>=1.0e-15) 
          regionIndex[number++] = i;
        else
          region[i]=0.0;
      }
      }
      regionSparse->setNumElements(number);
    }
  }
#endif
}
// Updates part of column (FTRANL) when densish
void 
CoinFactorization::updateColumnLDensish ( CoinIndexedVector * regionSparse ,
                                int * COIN_RESTRICT regionIndex)
  const
{
  double * COIN_RESTRICT region = regionSparse->denseVector (  );
  int number = regionSparse->getNumElements (  );
  int numberNonZero;
  double tolerance = zeroTolerance_;
  
  numberNonZero = 0;
  
  const CoinBigIndex * COIN_RESTRICT startColumn = startColumnL_.array();
  const int * COIN_RESTRICT indexRow = indexRowL_.array();
  const CoinFactorizationDouble * COIN_RESTRICT element = elementL_.array();
  int last = numberRows_;
  assert ( last == baseL_ + numberL_);
#if DENSE_CODE==1
  //can take out last bit of sparse L as empty
  last -= numberDense_;
#endif
  int smallestIndex = numberRowsExtra_;
  // do easy ones
  for (int k=0;k<number;k++) {
    int iPivot=regionIndex[k];
    if (iPivot>=baseL_) 
      smallestIndex = CoinMin(iPivot,smallestIndex);
    else
      regionIndex[numberNonZero++]=iPivot;
  }
  // now others
  for (int i = smallestIndex; i < last; i++ ) {
    CoinFactorizationDouble pivotValue = region[i];
    
    if ( fabs(pivotValue) > tolerance ) {
      CoinBigIndex start = startColumn[i];
      CoinBigIndex end = startColumn[i + 1];
      for (CoinBigIndex j = start; j < end; j ++ ) {
      int iRow = indexRow[j];
      CoinFactorizationDouble result = region[iRow];
      CoinFactorizationDouble value = element[j];

      region[iRow] = result - value * pivotValue;
      }     
      regionIndex[numberNonZero++] = i;
    } else {
      region[i] = 0.0;
    }       
  }     
  // and dense
  for (int i=last ; i < numberRows_; i++ ) {
    CoinFactorizationDouble pivotValue = region[i];
    if ( fabs(pivotValue) > tolerance ) {
      regionIndex[numberNonZero++] = i;
    } else {
      region[i] = 0.0;
    }       
  }     
  regionSparse->setNumElements ( numberNonZero );
} 
// Updates part of column (FTRANL) when sparsish
void 
CoinFactorization::updateColumnLSparsish ( CoinIndexedVector * regionSparse,
                                 int * COIN_RESTRICT regionIndex)
  const
{
  double * COIN_RESTRICT region = regionSparse->denseVector (  );
  int number = regionSparse->getNumElements (  );
  int numberNonZero;
  double tolerance = zeroTolerance_;
  
  numberNonZero = 0;
  
  const CoinBigIndex *startColumn = startColumnL_.array();
  const int *indexRow = indexRowL_.array();
  const CoinFactorizationDouble *element = elementL_.array();
  int last = numberRows_;
  assert ( last == baseL_ + numberL_);
#if DENSE_CODE==1
  //can take out last bit of sparse L as empty
  last -= numberDense_;
#endif
  // mark known to be zero
  int nInBig = sizeof(CoinBigIndex)/sizeof(int);
  CoinCheckZero * COIN_RESTRICT mark = reinterpret_cast<CoinCheckZero *> (sparse_.array()+(2+nInBig)*maximumRowsExtra_);
  int smallestIndex = numberRowsExtra_;
  // do easy ones
  for (int k=0;k<number;k++) {
    int iPivot=regionIndex[k];
    if (iPivot<baseL_) { 
      regionIndex[numberNonZero++]=iPivot;
    } else {
      smallestIndex = CoinMin(iPivot,smallestIndex);
      int iWord = iPivot>>CHECK_SHIFT;
      int iBit = iPivot-(iWord<<CHECK_SHIFT);
      if (mark[iWord]) {
      mark[iWord] = static_cast<CoinCheckZero>(mark[iWord] | (1<<iBit));
      } else {
      mark[iWord] = static_cast<CoinCheckZero>(1<<iBit);
      }
    }
  }
  // now others
  // First do up to convenient power of 2
  int jLast = (smallestIndex+BITS_PER_CHECK-1)>>CHECK_SHIFT;
  jLast = CoinMin((jLast<<CHECK_SHIFT),last);
  int i;
  for ( i = smallestIndex; i < jLast; i++ ) {
    CoinFactorizationDouble pivotValue = region[i];
    CoinBigIndex start = startColumn[i];
    CoinBigIndex end = startColumn[i + 1];
    
    if ( fabs(pivotValue) > tolerance ) {
      for (CoinBigIndex j = start; j < end; j ++ ) {
      int iRow = indexRow[j];
      CoinFactorizationDouble result = region[iRow];
      CoinFactorizationDouble value = element[j];
      region[iRow] = result - value * pivotValue;
      int iWord = iRow>>CHECK_SHIFT;
      int iBit = iRow-(iWord<<CHECK_SHIFT);
      if (mark[iWord]) {
        mark[iWord] = static_cast<CoinCheckZero>(mark[iWord] | (1<<iBit));
      } else {
        mark[iWord] = static_cast<CoinCheckZero>(1<<iBit);
      }
      }     
      regionIndex[numberNonZero++] = i;
    } else {
      region[i] = 0.0;
    }       
  }
  
  int kLast = last>>CHECK_SHIFT;
  if (jLast<last) {
    // now do in chunks
    for (int k=(jLast>>CHECK_SHIFT);k<kLast;k++) {
      unsigned int iMark = mark[k];
      if (iMark) {
      // something in chunk - do all (as imark may change)
      i = k<<CHECK_SHIFT;
      int iLast = i+BITS_PER_CHECK;
      for ( ; i < iLast; i++ ) {
        CoinFactorizationDouble pivotValue = region[i];
        CoinBigIndex start = startColumn[i];
        CoinBigIndex end = startColumn[i + 1];
        
        if ( fabs(pivotValue) > tolerance ) {
          CoinBigIndex j;
          for ( j = start; j < end; j ++ ) {
            int iRow = indexRow[j];
            CoinFactorizationDouble result = region[iRow];
            CoinFactorizationDouble value = element[j];
            region[iRow] = result - value * pivotValue;
            int iWord = iRow>>CHECK_SHIFT;
            int iBit = iRow-(iWord<<CHECK_SHIFT);
            if (mark[iWord]) {
            mark[iWord] = static_cast<CoinCheckZero>(mark[iWord] | (1<<iBit));
            } else {
            mark[iWord] = static_cast<CoinCheckZero>(1<<iBit);
            }
          }     
          regionIndex[numberNonZero++] = i;
        } else {
          region[i] = 0.0;
        }       
      }
      mark[k]=0; // zero out marked
      }
    }
    i = kLast<<CHECK_SHIFT;
  }
  for ( ; i < last; i++ ) {
    CoinFactorizationDouble pivotValue = region[i];
    CoinBigIndex start = startColumn[i];
    CoinBigIndex end = startColumn[i + 1];
    
    if ( fabs(pivotValue) > tolerance ) {
      for (CoinBigIndex j = start; j < end; j ++ ) {
      int iRow = indexRow[j];
      CoinFactorizationDouble result = region[iRow];
      CoinFactorizationDouble value = element[j];
      region[iRow] = result - value * pivotValue;
      }     
      regionIndex[numberNonZero++] = i;
    } else {
      region[i] = 0.0;
    }       
  }
  // Now dense part
  for ( ; i < numberRows_; i++ ) {
    double pivotValue = region[i];
    if ( fabs(pivotValue) > tolerance ) {
      regionIndex[numberNonZero++] = i;
    } else {
      region[i] = 0.0;
    }       
  }
  // zero out ones that might have been skipped
  mark[smallestIndex>>CHECK_SHIFT]=0;
  int kkLast = (numberRows_+BITS_PER_CHECK-1)>>CHECK_SHIFT;
  CoinZeroN(mark+kLast,kkLast-kLast);
  regionSparse->setNumElements ( numberNonZero );
}
// Updates part of column (FTRANL) when sparse
void 
CoinFactorization::updateColumnLSparse ( CoinIndexedVector * regionSparse ,
                                 int * COIN_RESTRICT regionIndex)
  const
{
  double * COIN_RESTRICT region = regionSparse->denseVector (  );
  int number = regionSparse->getNumElements (  );
  int numberNonZero;
  double tolerance = zeroTolerance_;
  
  numberNonZero = 0;
  
  const CoinBigIndex *startColumn = startColumnL_.array();
  const int *indexRow = indexRowL_.array();
  const CoinFactorizationDouble *element = elementL_.array();
  // use sparse_ as temporary area
  // mark known to be zero
  int * COIN_RESTRICT stack = sparse_.array();  /* pivot */
  int * COIN_RESTRICT list = stack + maximumRowsExtra_;  /* final list */
  CoinBigIndex * COIN_RESTRICT next = reinterpret_cast<CoinBigIndex *> (list + maximumRowsExtra_);  /* jnext */
  char * COIN_RESTRICT mark = reinterpret_cast<char *> (next + maximumRowsExtra_);
  int nList;
#ifdef COIN_DEBUG
  for (int i=0;i<maximumRowsExtra_;i++) {
    assert (!mark[i]);
  }
#endif
  nList=0;
  for (int k=0;k<number;k++) {
    int kPivot=regionIndex[k];
    if (kPivot>=baseL_) {
      assert (kPivot<numberRowsExtra_);
      //if (kPivot>=numberRowsExtra_) abort();
      if(!mark[kPivot]) {
      stack[0]=kPivot;
      CoinBigIndex j=startColumn[kPivot+1]-1;
        int nStack=0;
      while (nStack>=0) {
        /* take off stack */
        if (j>=startColumn[kPivot]) {
          int jPivot=indexRow[j--];
          assert (jPivot>=baseL_&&jPivot<numberRowsExtra_);
          //if (jPivot<baseL_||jPivot>=numberRowsExtra_) abort();
          /* put back on stack */
          next[nStack] =j;
          if (!mark[jPivot]) {
            /* and new one */
            kPivot=jPivot;
            j = startColumn[kPivot+1]-1;
            stack[++nStack]=kPivot;
            assert (kPivot<numberRowsExtra_);
            //if (kPivot>=numberRowsExtra_) abort();
            mark[kPivot]=1;
            next[nStack]=j;
          }
        } else {
          /* finished so mark */
          list[nList++]=kPivot;
          mark[kPivot]=1;
          --nStack;
          if (nStack>=0) {
            kPivot=stack[nStack];
            assert (kPivot<numberRowsExtra_);
            j=next[nStack];
          }
        }
      }
      }
    } else {
      // just put on list
      regionIndex[numberNonZero++]=kPivot;
    }
  }
  for (int i=nList-1;i>=0;i--) {
    int iPivot = list[i];
    mark[iPivot]=0;
    CoinFactorizationDouble pivotValue = region[iPivot];
    if ( fabs ( pivotValue ) > tolerance ) {
      regionIndex[numberNonZero++]=iPivot;
      for (CoinBigIndex j = startColumn[iPivot]; 
         j < startColumn[iPivot+1]; j ++ ) {
      int iRow = indexRow[j];
      CoinFactorizationDouble value = element[j];
      region[iRow] -= value * pivotValue;
      }
    } else {
      region[iPivot]=0.0;
    }
  }
  regionSparse->setNumElements ( numberNonZero );
}
/* Updates one column (FTRAN) from region2
   Tries to do FT update
   number returned is negative if no room.
   Also updates region3
   region1 starts as zero and is zero at end */
int 
00537 CoinFactorization::updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
                              CoinIndexedVector * regionSparse2,
                              CoinIndexedVector * regionSparse3,
                              bool noPermuteRegion3)
{
#if 1
  //#ifdef NDEBUG
  //#undef NDEBUG
  //#endif
  //#define COIN_DEBUG
#ifdef COIN_DEBUG
  regionSparse1->checkClean();
  CoinIndexedVector save2(*regionSparse2);
  CoinIndexedVector save3(*regionSparse3);
#endif
  CoinIndexedVector * regionFT ;
  CoinIndexedVector * regionUpdate ;
  int * COIN_RESTRICT regionIndex ;
  int numberNonZero ;
  const int *permute = permute_.array();
  int * COIN_RESTRICT index ;
  double * COIN_RESTRICT region ;
  if (!noPermuteRegion3) {
    regionFT = regionSparse3;
    regionUpdate = regionSparse1;
    //permute and move indices into index array
    regionIndex = regionUpdate->getIndices (  );
    //int numberNonZero;
    region = regionUpdate->denseVector();
    
    numberNonZero = regionSparse3->getNumElements();
    int * COIN_RESTRICT index = regionSparse3->getIndices();
    double * COIN_RESTRICT array = regionSparse3->denseVector();
    assert (!regionSparse3->packedMode());
    for (int j = 0; j < numberNonZero; j ++ ) {
      int iRow = index[j];
      double value = array[iRow];
      array[iRow]=0.0;
      iRow = permute[iRow];
      region[iRow] = value;
      regionIndex[j] = iRow;
    }
    regionUpdate->setNumElements ( numberNonZero );
  } else {
    regionFT = regionSparse1;
    regionUpdate = regionSparse3;
  }
  //permute and move indices into index array (in U)
  regionIndex = regionFT->getIndices (  );
  numberNonZero = regionSparse2->getNumElements();
  index = regionSparse2->getIndices();
  region = regionFT->denseVector();
  double * COIN_RESTRICT array = regionSparse2->denseVector();
  CoinBigIndex * COIN_RESTRICT startColumnU = startColumnU_.array();
  CoinBigIndex start = startColumnU[maximumColumnsExtra_];
  startColumnU[numberColumnsExtra_] = start;
  regionIndex = indexRowU_.array() + start;

  assert(regionSparse2->packedMode());
  for (int j = 0; j < numberNonZero; j ++ ) {
    int iRow = index[j];
    double value = array[j];
    array[j]=0.0;
    iRow = permute[iRow];
    region[iRow] = value;
    regionIndex[j] = iRow;
  }
  regionFT->setNumElements ( numberNonZero );
  if (collectStatistics_) {
    numberFtranCounts_+=2;
    ftranCountInput_ += regionFT->getNumElements()+
      regionUpdate->getNumElements();
  }
    
  //  ******* L
  updateColumnL ( regionFT, regionIndex );
  updateColumnL ( regionUpdate, regionUpdate->getIndices() );
  if (collectStatistics_) 
    ftranCountAfterL_ += regionFT->getNumElements()+
      regionUpdate->getNumElements();
  //permute extra
  //row bits here
  updateColumnRFT ( regionFT, regionIndex );
  updateColumnR ( regionUpdate );
  if (collectStatistics_) 
    ftranCountAfterR_ += regionFT->getNumElements()+
    regionUpdate->getNumElements();
  //  ******* U - see if densish
  // Guess at number at end
  int goSparse=0;
  if (sparseThreshold_>0) {
    int numberNonZero = (regionUpdate->getNumElements (  )+
                   regionFT->getNumElements())>>1;
    if (ftranAverageAfterR_) {
      int newNumber = static_cast<int> (numberNonZero*ftranAverageAfterU_);
      if (newNumber< sparseThreshold_)
      goSparse = 2;
      else if (newNumber< sparseThreshold2_)
      goSparse = 1;
    } else {
      if (numberNonZero<sparseThreshold_) 
      goSparse = 2;
    }
  }
#ifndef COIN_FAST_CODE
  assert (slackValue_==-1.0);
#endif
  if (!goSparse&&numberRows_<1000) {
    double * COIN_RESTRICT arrayFT = regionFT->denseVector (  );
    int * COIN_RESTRICT indexFT = regionFT->getIndices();
    int numberNonZeroFT;
    double * COIN_RESTRICT arrayUpdate = regionUpdate->denseVector (  );
    int * COIN_RESTRICT indexUpdate = regionUpdate->getIndices();
    int numberNonZeroUpdate;
    updateTwoColumnsUDensish(numberNonZeroFT,arrayFT,indexFT,
                       numberNonZeroUpdate,arrayUpdate,indexUpdate);
    regionFT->setNumElements ( numberNonZeroFT );
    regionUpdate->setNumElements ( numberNonZeroUpdate );
  } else {
    // sparse 
    updateColumnU ( regionFT, regionIndex);
    updateColumnU ( regionUpdate, regionUpdate->getIndices());
  }
  permuteBack(regionFT,regionSparse2);
  if (!noPermuteRegion3) {
    permuteBack(regionUpdate,regionSparse3);
  }
#ifdef COIN_DEBUG
  int n2=regionSparse2->getNumElements();
  regionSparse1->checkClean();
  int n2a= updateColumnFT(regionSparse1,&save2);
  assert(n2==n2a);
  {
    int j;
    double * regionA = save2.denseVector();
    int * indexA = save2.getIndices();
    double * regionB = regionSparse2->denseVector();
    int * indexB = regionSparse2->getIndices();
    for (j=0;j<n2;j++) {
      int k = indexA[j];
      assert (k==indexB[j]);
      CoinFactorizationDouble value = regionA[j];
      assert (value==regionB[j]);
    }
  }
  updateColumn(&save3,
             &save3,
             noPermuteRegion3);
  int n3=regionSparse3->getNumElements();
  assert (n3==save3.getNumElements());
  {
    int j;
    double * regionA = save3.denseVector();
    int * indexA = save3.getIndices();
    double * regionB = regionSparse3->denseVector();
    int * indexB = regionSparse3->getIndices();
    for (j=0;j<n3;j++) {
      int k = indexA[j];
      assert (k==indexB[j]);
      CoinFactorizationDouble value = regionA[k];
      assert (value==regionB[k]);
    }
  }
  //*regionSparse2=save2;
  //*regionSparse3=save3;
  printf("REGION2 %d els\n",regionSparse2->getNumElements());
  regionSparse2->print();
  printf("REGION3 %d els\n",regionSparse3->getNumElements());
  regionSparse3->print();
#endif
  return regionSparse2->getNumElements();
#else
  int returnCode= updateColumnFT(regionSparse1,
                         regionSparse2);
  assert (noPermuteRegion3);
  updateColumn(regionSparse3,
             regionSparse3,
             noPermuteRegion3);
  //printf("REGION2 %d els\n",regionSparse2->getNumElements());
  //regionSparse2->print();
  //printf("REGION3 %d els\n",regionSparse3->getNumElements());
  //regionSparse3->print();
  return returnCode;
#endif
}
// Updates part of 2 columns (FTRANU) real work
void
00724 CoinFactorization::updateTwoColumnsUDensish (
                                   int & numberNonZero1,
                                   double * COIN_RESTRICT region1, 
                                   int * COIN_RESTRICT index1,
                                   int & numberNonZero2,
                                   double * COIN_RESTRICT region2, 
                                   int * COIN_RESTRICT index2) const
{
  double tolerance = zeroTolerance_;
  const CoinBigIndex * COIN_RESTRICT startColumn = startColumnU_.array();
  const int * COIN_RESTRICT indexRow = indexRowU_.array();
  const CoinFactorizationDouble * COIN_RESTRICT element = elementU_.array();
  int numberNonZeroA = 0;
  int numberNonZeroB = 0;
  const int *numberInColumn = numberInColumn_.array();
  const CoinFactorizationDouble *pivotRegion = pivotRegion_.array();
  
  for (int i = numberU_-1 ; i >= numberSlacks_; i-- ) {
    CoinFactorizationDouble pivotValue2 = region2[i];
    region2[i] = 0.0;
    CoinFactorizationDouble pivotValue1 = region1[i];
    region1[i] = 0.0;
    if ( fabs ( pivotValue2 ) > tolerance ) {
      CoinBigIndex start = startColumn[i];
      const CoinFactorizationDouble * COIN_RESTRICT thisElement = element+start;
      const int * COIN_RESTRICT thisIndex = indexRow+start;
      if ( fabs ( pivotValue1 ) <= tolerance ) {
      // just region 2
      for (CoinBigIndex j=numberInColumn[i]-1 ; j >=0; j-- ) {
        int iRow = thisIndex[j];
        CoinFactorizationDouble value = thisElement[j];
#ifdef NO_LOAD
        region2[iRow] -= value * pivotValue2;
#else
        CoinFactorizationDouble regionValue2 = region2[iRow];
        region2[iRow] = regionValue2 - value * pivotValue2;
#endif
      }
      pivotValue2 *= pivotRegion[i];
      region2[i]=pivotValue2;
      index2[numberNonZeroB++]=i;
      } else {
      // both
      for (CoinBigIndex j=numberInColumn[i]-1 ; j >=0; j-- ) {
        int iRow = thisIndex[j];
        CoinFactorizationDouble value = thisElement[j];
#ifdef NO_LOAD
        region1[iRow] -= value * pivotValue1;
        region2[iRow] -= value * pivotValue2;
#else
        CoinFactorizationDouble regionValue1 = region1[iRow];
        CoinFactorizationDouble regionValue2 = region2[iRow];
        region1[iRow] = regionValue1 - value * pivotValue1;
        region2[iRow] = regionValue2 - value * pivotValue2;
#endif
      }
      pivotValue1 *= pivotRegion[i];
      pivotValue2 *= pivotRegion[i];
      region1[i]=pivotValue1;
      index1[numberNonZeroA++]=i;
      region2[i]=pivotValue2;
      index2[numberNonZeroB++]=i;
      }
    } else if ( fabs ( pivotValue1 ) > tolerance ) {
      CoinBigIndex start = startColumn[i];
      const CoinFactorizationDouble * COIN_RESTRICT thisElement = element+start;
      const int * COIN_RESTRICT thisIndex = indexRow+start;
      // just region 1
      for (CoinBigIndex j=numberInColumn[i]-1 ; j >=0; j-- ) {
      int iRow = thisIndex[j];
      CoinFactorizationDouble value = thisElement[j];
#ifdef NO_LOAD
      region1[iRow] -= value * pivotValue1;
#else
      CoinFactorizationDouble regionValue1 = region1[iRow];
      region1[iRow] = regionValue1 - value * pivotValue1;
#endif
      }
      pivotValue1 *= pivotRegion[i];
      region1[i]=pivotValue1;
      index1[numberNonZeroA++]=i;
    }
  }
  // Slacks 
    
  for (int i = numberSlacks_-1; i>=0;i--) {
    double value2 = region2[i];
    double value1 = region1[i];
    bool value1NonZero = (value1!=0.0);
    if ( fabs(value2) > tolerance ) {
      region2[i]=-value2;
      index2[numberNonZeroB++]=i;
    } else {
      region2[i]=0.0;
    }
    if ( value1NonZero ) {
      index1[numberNonZeroA]=i;
      if ( fabs(value1) > tolerance ) {
      region1[i]=-value1;
      numberNonZeroA++;
      } else {
      region1[i]=0.0;
      }
    }
  }
  numberNonZero1=numberNonZeroA;
  numberNonZero2=numberNonZeroB;
}
//  updateColumnU.  Updates part of column (FTRANU)
void
00834 CoinFactorization::updateColumnU ( CoinIndexedVector * regionSparse,
                           int * indexIn) const
{
  int numberNonZero = regionSparse->getNumElements (  );

  int goSparse;
  // Guess at number at end
  if (sparseThreshold_>0) {
    if (ftranAverageAfterR_) {
      int newNumber = static_cast<int> (numberNonZero*ftranAverageAfterU_);
      if (newNumber< sparseThreshold_)
      goSparse = 2;
      else if (newNumber< sparseThreshold2_)
      goSparse = 1;
      else
      goSparse = 0;
    } else {
      if (numberNonZero<sparseThreshold_) 
      goSparse = 2;
      else
      goSparse = 0;
    }
  } else {
    goSparse=0;
  }
  switch (goSparse) {
  case 0: // densish
    {
      double *region = regionSparse->denseVector (  );
      int * regionIndex = regionSparse->getIndices();
      int numberNonZero=updateColumnUDensish(region,regionIndex);
      regionSparse->setNumElements ( numberNonZero );
    }
    break;
  case 1: // middling
    updateColumnUSparsish(regionSparse,indexIn);
    break;
  case 2: // sparse
    updateColumnUSparse(regionSparse,indexIn);
    break;
  }
  if (collectStatistics_) 
    ftranCountAfterU_ += regionSparse->getNumElements (  );
}
#ifdef COIN_DEVELOP
double ncall_DZ=0.0;
double nrow_DZ=0.0;
double nslack_DZ=0.0;
double nU_DZ=0.0;
double nnz_DZ=0.0;
double nDone_DZ=0.0;
#endif
// Updates part of column (FTRANU) real work
int 
00888 CoinFactorization::updateColumnUDensish ( double * COIN_RESTRICT region, 
                                int * COIN_RESTRICT regionIndex) const
{
  double tolerance = zeroTolerance_;
  const CoinBigIndex *startColumn = startColumnU_.array();
  const int *indexRow = indexRowU_.array();
  const CoinFactorizationDouble *element = elementU_.array();
  int numberNonZero = 0;
  const int *numberInColumn = numberInColumn_.array();
  const CoinFactorizationDouble *pivotRegion = pivotRegion_.array();
#ifdef COIN_DEVELOP
  ncall_DZ++;
  nrow_DZ += numberRows_;
  nslack_DZ += numberSlacks_;
  nU_DZ += numberU_;
#endif
  
  for (int i = numberU_-1 ; i >= numberSlacks_; i-- ) {
    CoinFactorizationDouble pivotValue = region[i];
    if (pivotValue) {
#ifdef COIN_DEVELOP
      nnz_DZ++;
#endif
      region[i] = 0.0;
      if ( fabs ( pivotValue ) > tolerance ) {
      CoinBigIndex start = startColumn[i];
      const CoinFactorizationDouble * thisElement = element+start;
      const int * thisIndex = indexRow+start;
#ifdef COIN_DEVELOP
      nDone_DZ += numberInColumn[i];
#endif
      for (CoinBigIndex j=numberInColumn[i]-1 ; j >=0; j-- ) {
        int iRow = thisIndex[j];
        CoinFactorizationDouble regionValue = region[iRow];
        CoinFactorizationDouble value = thisElement[j];
        region[iRow] = regionValue - value * pivotValue;
      }
      pivotValue *= pivotRegion[i];
      region[i]=pivotValue;
      regionIndex[numberNonZero++]=i;
      }
    }
  }
    
  // now do slacks
#ifndef COIN_FAST_CODE
  if (slackValue_==-1.0) {
#endif
#if 0
    // Could skew loop to pick up next one earlier
    // might improve pipelining
    for (int i = numberSlacks_-1; i>2;i-=2) {
      double value0 = region[i];
      double absValue0 = fabs ( value0 );
      double value1 = region[i-1];
      double absValue1 = fabs ( value1 );
      if ( value0 ) {
      if ( absValue0 > tolerance ) {
        region[i]=-value0;
        regionIndex[numberNonZero++]=i;
      } else {
        region[i]=0.0;
      }
      }
      if ( value1 ) {
      if ( absValue1 > tolerance ) {
        region[i-1]=-value1;
        regionIndex[numberNonZero++]=i-1;
      } else {
        region[i-1]=0.0;
      }
      }
    }
    for ( ; i>=0;i--) {
      double value = region[i];
      double absValue = fabs ( value );
      if ( value ) {
      if ( absValue > tolerance ) {
        region[i]=-value;
        regionIndex[numberNonZero++]=i;
      } else {
        region[i]=0.0;
      }
      }
    }
#else
    for (int i = numberSlacks_-1; i>=0;i--) {
      double value = region[i];
      if ( value ) {
      region[i]=-value;
      regionIndex[numberNonZero]=i;
      if ( fabs(value) > tolerance ) 
        numberNonZero++;
      else 
        region[i]=0.0;
      }
    }
#endif
#ifndef COIN_FAST_CODE
  } else {
    assert (slackValue_==1.0);
    for (int i = numberSlacks_-1; i>=0;i--) {
      double value = region[i];
      double absValue = fabs ( value );
      if ( value ) {
      region[i]=0.0;
      if ( absValue > tolerance ) {
        region[i]=value;
        regionIndex[numberNonZero++]=i;
      }
      }
    }
  }
#endif
  return numberNonZero;
}
//  updateColumnU.  Updates part of column (FTRANU)
/*
  Since everything is in order I should be able to do a better job of
  marking stuff - think.  Also as L is static maybe I can do something
  better there (I know I could if I marked the depth of every element
  but that would lead to other inefficiencies.
*/
void
CoinFactorization::updateColumnUSparse ( CoinIndexedVector * regionSparse,
                               int * COIN_RESTRICT indexIn) const
{
  int numberNonZero = regionSparse->getNumElements (  );
  int * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
  double * COIN_RESTRICT region = regionSparse->denseVector (  );
  double tolerance = zeroTolerance_;
  const CoinBigIndex *startColumn = startColumnU_.array();
  const int *indexRow = indexRowU_.array();
  const CoinFactorizationDouble *element = elementU_.array();
  const CoinFactorizationDouble *pivotRegion = pivotRegion_.array();
  // use sparse_ as temporary area
  // mark known to be zero
  int * COIN_RESTRICT stack = sparse_.array();  /* pivot */
  int * COIN_RESTRICT list = stack + maximumRowsExtra_;  /* final list */
  CoinBigIndex * COIN_RESTRICT next = reinterpret_cast<CoinBigIndex *> (list + maximumRowsExtra_);  /* jnext */
  char * COIN_RESTRICT mark = reinterpret_cast<char *> (next + maximumRowsExtra_);
#ifdef COIN_DEBUG
  for (int i=0;i<maximumRowsExtra_;i++) {
    assert (!mark[i]);
  }
#endif

  // move slacks to end of stack list
  int * COIN_RESTRICT putLast = stack+maximumRowsExtra_;
  int * COIN_RESTRICT put = putLast;

  const int *numberInColumn = numberInColumn_.array();
  int nList = 0;
  for (int i=0;i<numberNonZero;i++) {
    int kPivot=indexIn[i];
    stack[0]=kPivot;
    CoinBigIndex j=startColumn[kPivot]+numberInColumn[kPivot]-1;
    int nStack=1;
    next[0]=j;
    while (nStack) {
      /* take off stack */
      int kPivot=stack[--nStack];
      if (mark[kPivot]!=1) {
      j=next[nStack];
      if (j>=startColumn[kPivot]) {
        kPivot=indexRow[j--];
        /* put back on stack */
        next[nStack++] =j;
        if (!mark[kPivot]) {
          /* and new one */
          int numberIn = numberInColumn[kPivot];
          if (numberIn) {
            j = startColumn[kPivot]+numberIn-1;
            stack[nStack]=kPivot;
            mark[kPivot]=2;
            next[nStack++]=j;
          } else {
            // can do immediately
            /* finished so mark */
            mark[kPivot]=1;
            if (kPivot>=numberSlacks_) {
            list[nList++]=kPivot;
            } else {
            // slack - put at end
            --put;
            *put=kPivot;
            }
          }
        }
      } else {
        /* finished so mark */
        mark[kPivot]=1;
        if (kPivot>=numberSlacks_) {
          list[nList++]=kPivot;
        } else {
          // slack - put at end
          assert (!numberInColumn[kPivot]);
          --put;
          *put=kPivot;
        }
      }
      }
    }
  }
#if 0
  {
    std::sort(list,list+nList);
    int i;
    int last;
    last =-1;
    for (i=0;i<nList;i++) {
      int k = list[i];
      assert (k>last);
      last=k;
    }
    std::sort(put,putLast);
    int n = putLast-put;
    last =-1;
    for (i=0;i<n;i++) {
      int k = put[i];
      assert (k>last);
      last=k;
    }
  }
#endif
  numberNonZero=0;
  for (int i=nList-1;i>=0;i--) {
    int iPivot = list[i];
    mark[iPivot]=0;
    CoinFactorizationDouble pivotValue = region[iPivot];
    region[iPivot]=0.0;
    if ( fabs ( pivotValue ) > tolerance ) {
      CoinBigIndex start = startColumn[iPivot];
      int number = numberInColumn[iPivot];
      
      CoinBigIndex j;
      for ( j = start; j < start+number; j++ ) {
      CoinFactorizationDouble value = element[j];
      int iRow = indexRow[j];
      region[iRow] -=  value * pivotValue;
      }
      pivotValue *= pivotRegion[iPivot];
      region[iPivot]=pivotValue;
      regionIndex[numberNonZero++]=iPivot;
    }
  }
  // slacks
#ifndef COIN_FAST_CODE
  if (slackValue_==1.0) {
    for (;put<putLast;put++) {
      int iPivot = *put;
      mark[iPivot]=0;
      CoinFactorizationDouble pivotValue = region[iPivot];
      region[iPivot]=0.0;
      if ( fabs ( pivotValue ) > tolerance ) {
      region[iPivot]=pivotValue;
      regionIndex[numberNonZero++]=iPivot;
      }
    }
  } else {
#endif
    for (;put<putLast;put++) {
      int iPivot = *put;
      mark[iPivot]=0;
      CoinFactorizationDouble pivotValue = region[iPivot];
      region[iPivot]=0.0;
      if ( fabs ( pivotValue ) > tolerance ) {
      region[iPivot]=-pivotValue;
      regionIndex[numberNonZero++]=iPivot;
      }
    }
#ifndef COIN_FAST_CODE
  }
#endif
  regionSparse->setNumElements ( numberNonZero );
}
//  updateColumnU.  Updates part of column (FTRANU)
/*
  Since everything is in order I should be able to do a better job of
  marking stuff - think.  Also as L is static maybe I can do something
  better there (I know I could if I marked the depth of every element
  but that would lead to other inefficiencies.
*/
#ifdef COIN_DEVELOP
double ncall_SZ=0.0;
double nrow_SZ=0.0;
double nslack_SZ=0.0;
double nU_SZ=0.0;
double nnz_SZ=0.0;
double nDone_SZ=0.0;
#endif
void
CoinFactorization::updateColumnUSparsish ( CoinIndexedVector * regionSparse,
                                 int * COIN_RESTRICT indexIn) const
{
  int * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
  // mark known to be zero
  int * COIN_RESTRICT stack = sparse_.array();  /* pivot */
  int * COIN_RESTRICT list = stack + maximumRowsExtra_;  /* final list */
  CoinBigIndex * COIN_RESTRICT next = reinterpret_cast<CoinBigIndex *> (list + maximumRowsExtra_);  /* jnext */
  CoinCheckZero * COIN_RESTRICT mark = reinterpret_cast<CoinCheckZero *> (next + maximumRowsExtra_);
  const int *numberInColumn = numberInColumn_.array();
#ifdef COIN_DEBUG
  for (int i=0;i<maximumRowsExtra_;i++) {
    assert (!mark[i]);
  }
#endif

  int nMarked=0;
  int numberNonZero = regionSparse->getNumElements (  );
  double * COIN_RESTRICT region = regionSparse->denseVector (  );
  double tolerance = zeroTolerance_;
  const CoinBigIndex *startColumn = startColumnU_.array();
  const int *indexRow = indexRowU_.array();
  const CoinFactorizationDouble *element = elementU_.array();
  const CoinFactorizationDouble *pivotRegion = pivotRegion_.array();
#ifdef COIN_DEVELOP
  ncall_SZ++;
  nrow_SZ += numberRows_;
  nslack_SZ += numberSlacks_;
  nU_SZ += numberU_;
#endif

  // Bjarni: renamed i to ii to avoid in MSVC compiler potential conflicts warnings with i defined in line 1230
  for (int ii = 0; ii < numberNonZero; ii++) { 
    int iPivot=indexIn[ii];
    int iWord = iPivot>>CHECK_SHIFT;
    int iBit = iPivot-(iWord<<CHECK_SHIFT);
    if (mark[iWord]) {
      mark[iWord] = static_cast<CoinCheckZero>(mark[iWord] | (1<<iBit));
    } else {
      mark[iWord] = static_cast<CoinCheckZero>(1<<iBit);
      stack[nMarked++]=iWord;
    }
  }
  numberNonZero = 0;
  // First do down to convenient power of 2
  CoinBigIndex jLast = (numberU_-1)>>CHECK_SHIFT;
  jLast = CoinMax((jLast<<CHECK_SHIFT),static_cast<CoinBigIndex> (numberSlacks_));
  int i;
  for ( i = numberU_-1 ; i >= jLast; i-- ) {
    CoinFactorizationDouble pivotValue = region[i];
    region[i] = 0.0;
    if ( fabs ( pivotValue ) > tolerance ) {
#ifdef COIN_DEVELOP
      nnz_SZ ++;
#endif
      CoinBigIndex start = startColumn[i];
      const CoinFactorizationDouble * thisElement = element+start;
      const int * thisIndex = indexRow+start;
      
#ifdef COIN_DEVELOP
      nDone_SZ += numberInColumn[i];
#endif
      for (int j=numberInColumn[i]-1 ; j >=0; j-- ) {
      int iRow0 = thisIndex[j];
      CoinFactorizationDouble regionValue0 = region[iRow0];
      CoinFactorizationDouble value0 = thisElement[j];
      int iWord = iRow0>>CHECK_SHIFT;
      int iBit = iRow0-(iWord<<CHECK_SHIFT);
      if (mark[iWord]) {
        mark[iWord] = static_cast<CoinCheckZero>(mark[iWord] | (1<<iBit));
      } else {
        mark[iWord] = static_cast<CoinCheckZero>(1<<iBit);
        stack[nMarked++]=iWord;
      }
      region[iRow0] = regionValue0 - value0 * pivotValue;
      }
      pivotValue *= pivotRegion[i];
      region[i]=pivotValue;
      regionIndex[numberNonZero++]=i;
    }
  }
  int kLast = (numberSlacks_+BITS_PER_CHECK-1)>>CHECK_SHIFT;
  if (jLast>numberSlacks_) {
    // now do in chunks
    for (int k=(jLast>>CHECK_SHIFT)-1;k>=kLast;k--) {
      unsigned int iMark = mark[k];
      if (iMark) {
      // something in chunk - do all (as imark may change)
      int iLast = k<<CHECK_SHIFT;
      for ( i = iLast+BITS_PER_CHECK-1 ; i >= iLast; i-- ) {
        CoinFactorizationDouble pivotValue = region[i];
        if (pivotValue) {
#ifdef COIN_DEVELOP
          nnz_SZ ++;
#endif
          region[i] = 0.0;
          if ( fabs ( pivotValue ) > tolerance ) {
            CoinBigIndex start = startColumn[i];
            const CoinFactorizationDouble * thisElement = element+start;
            const int * thisIndex = indexRow+start;
#ifdef COIN_DEVELOP
            nDone_SZ += numberInColumn[i];
#endif
            for (int j=numberInColumn[i]-1 ; j >=0; j-- ) {
            int iRow0 = thisIndex[j];
            CoinFactorizationDouble regionValue0 = region[iRow0];
            CoinFactorizationDouble value0 = thisElement[j];
            int iWord = iRow0>>CHECK_SHIFT;
            int iBit = iRow0-(iWord<<CHECK_SHIFT);
            if (mark[iWord]) {
              mark[iWord] = static_cast<CoinCheckZero>(mark[iWord] | (1<<iBit));
            } else {
              mark[iWord] = static_cast<CoinCheckZero>(1<<iBit);
              stack[nMarked++]=iWord;
            }
            region[iRow0] = regionValue0 - value0 * pivotValue;
            }
            pivotValue *= pivotRegion[i];
            region[i]=pivotValue;
            regionIndex[numberNonZero++]=i;
          }
        }
      }
      mark[k]=0;
      }
    }
    i = (kLast<<CHECK_SHIFT)-1;
  }
  for ( ; i >= numberSlacks_; i-- ) {
    CoinFactorizationDouble pivotValue = region[i];
    region[i] = 0.0;
    if ( fabs ( pivotValue ) > tolerance ) {
#ifdef COIN_DEVELOP
      nnz_SZ ++;
#endif
      CoinBigIndex start = startColumn[i];
      const CoinFactorizationDouble * thisElement = element+start;
      const int * thisIndex = indexRow+start;
#ifdef COIN_DEVELOP
      nDone_SZ += numberInColumn[i];
#endif
      for (int j=numberInColumn[i]-1 ; j >=0; j-- ) {
      int iRow0 = thisIndex[j];
      CoinFactorizationDouble regionValue0 = region[iRow0];
      CoinFactorizationDouble value0 = thisElement[j];
      int iWord = iRow0>>CHECK_SHIFT;
      int iBit = iRow0-(iWord<<CHECK_SHIFT);
      if (mark[iWord]) {
        mark[iWord] = static_cast<CoinCheckZero>(mark[iWord] | (1<<iBit));
      } else {
        mark[iWord] = static_cast<CoinCheckZero>(1<<iBit);
        stack[nMarked++]=iWord;
      }
      region[iRow0] = regionValue0 - value0 * pivotValue;
      }
      pivotValue *= pivotRegion[i];
      region[i]=pivotValue;
      regionIndex[numberNonZero++]=i;
    }
  }
  
  if (numberSlacks_) {
    // now do slacks
#ifndef COIN_FAST_CODE
    double factor = slackValue_;
    if (factor==1.0) {
      // First do down to convenient power of 2
      CoinBigIndex jLast = (numberSlacks_-1)>>CHECK_SHIFT;
      jLast = jLast<<CHECK_SHIFT;
      for ( i = numberSlacks_-1; i>=jLast;i--) {
      double value = region[i];
      double absValue = fabs ( value );
      if ( value ) {
        region[i]=0.0;
        if ( absValue > tolerance ) {
        region[i]=value;
        regionIndex[numberNonZero++]=i;
        }
      }
      }
      mark[jLast]=0;
      // now do in chunks
      for (int k=(jLast>>CHECK_SHIFT)-1;k>=0;k--) {
      unsigned int iMark = mark[k];
      if (iMark) {
        // something in chunk - do all (as imark may change)
        int iLast = k<<CHECK_SHIFT;
        i = iLast+BITS_PER_CHECK-1;
        for ( ; i >= iLast; i-- ) {
          double value = region[i];
          double absValue = fabs ( value );
          if ( value ) {
            region[i]=0.0;
            if ( absValue > tolerance ) {
            region[i]=value;
            regionIndex[numberNonZero++]=i;
            }
          }
        }
        mark[k]=0;
      }
      }
    } else {
      assert (factor==-1.0);
#endif
      // First do down to convenient power of 2
      CoinBigIndex jLast = (numberSlacks_-1)>>CHECK_SHIFT;
      jLast = jLast<<CHECK_SHIFT;
      for ( i = numberSlacks_-1; i>=jLast;i--) {
      double value = region[i];
      double absValue = fabs ( value );
      if ( value ) {
        region[i]=0.0;
        if ( absValue > tolerance ) {
          region[i]=-value;
          regionIndex[numberNonZero++]=i;
        }
      }
      }
      mark[jLast]=0;
      // now do in chunks
      for (int k=(jLast>>CHECK_SHIFT)-1;k>=0;k--) {
      unsigned int iMark = mark[k];
      if (iMark) {
        // something in chunk - do all (as imark may change)
        int iLast = k<<CHECK_SHIFT;
        i = iLast+BITS_PER_CHECK-1;
        for ( ; i >= iLast; i-- ) {
          double value = region[i];
          double absValue = fabs ( value );
          if ( value ) {
            region[i]=0.0;
            if ( absValue > tolerance ) {
            region[i]=-value;
            regionIndex[numberNonZero++]=i;
            }
          }
        }
        mark[k]=0;
      }
      }
#ifndef COIN_FAST_CODE
    }
#endif
  }
  regionSparse->setNumElements ( numberNonZero );
  mark[(numberU_-1)>>CHECK_SHIFT]=0;
  mark[numberSlacks_>>CHECK_SHIFT]=0;
  if (numberSlacks_)
    mark[(numberSlacks_-1)>>CHECK_SHIFT]=0;
#ifdef COIN_DEBUG
  for (i=0;i<maximumRowsExtra_;i++) {
    assert (!mark[i]);
  }
#endif
}
//  updateColumnR.  Updates part of column (FTRANR)
void
01438 CoinFactorization::updateColumnR ( CoinIndexedVector * regionSparse ) const
{
  double * COIN_RESTRICT region = regionSparse->denseVector (  );
  int * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
  int numberNonZero = regionSparse->getNumElements (  );

  if ( !numberR_ )
    return; //return if nothing to do
  double tolerance = zeroTolerance_;

  const CoinBigIndex * startColumn = startColumnR_.array()-numberRows_;
  const int * indexRow = indexRowR_;
  const CoinFactorizationDouble * element = elementR_;
  const int * permute = permute_.array();

  // Work out very dubious idea of what would be fastest
  int method=-1;
  // Size of R
  double sizeR=startColumnR_.array()[numberR_];
  // Average
  double averageR = sizeR/(static_cast<double> (numberRowsExtra_));
  // weights (relative to actual work)
  double setMark = 0.1; // setting mark
  double test1= 1.0; // starting ftran (without testPivot)
  double testPivot = 2.0; // Seeing if zero etc
  double startDot=2.0; // For starting dot product version
  // For final scan
  double final = numberNonZero*1.0;
  double methodTime[3];
  // For second type
  methodTime[1] = numberPivots_ * (testPivot + ((static_cast<double> (numberNonZero))/(static_cast<double> (numberRows_))
                                    * averageR));
  methodTime[1] += numberNonZero *(test1 + averageR);
  // For first type
  methodTime[0] = methodTime[1] + (numberNonZero+numberPivots_)*setMark;
  methodTime[1] += numberNonZero*final;
  // third
  methodTime[2] = sizeR + numberPivots_*startDot + numberNonZero*final;
  // switch off if necessary
  if (!numberInColumnPlus_.array()) {
    methodTime[0]=1.0e100;
    methodTime[1]=1.0e100;
  } else if (!sparse_.array()) {
    methodTime[0]=1.0e100;
  }
  double best=1.0e100;
  for (int i=0;i<3;i++) {
    if (methodTime[i]<best) {
      best=methodTime[i];
      method=i;
    }
  }
  assert (method>=0);
  const int * numberInColumnPlus = numberInColumnPlus_.array();
  //if (method==1)
  //printf(" methods %g %g %g - chosen %d\n",methodTime[0],methodTime[1],methodTime[2],method);

  switch (method) {
  case 0:
#ifdef STACK
    {
      // use sparse_ as temporary area
      // mark known to be zero
      int * COIN_RESTRICT stack = sparse_.array();  /* pivot */
      int * COIN_RESTRICT list = stack + maximumRowsExtra_;  /* final list */
      CoinBigIndex * COIN_RESTRICT next = (CoinBigIndex *) (list + maximumRowsExtra_);  /* jnext */
      char * COIN_RESTRICT mark = (char *) (next + maximumRowsExtra_);
      // we have another copy of R in R
      const CoinFactorizationDouble * elementR = elementR_ + lengthAreaR_;
      const int * indexRowR = indexRowR_ + lengthAreaR_;
      const CoinBigIndex * startR = startColumnR_.array()+maximumPivots_+1;
      int nList=0;
      const int * permuteBack = permuteBack_.array();
      for (int k=0;k<numberNonZero;k++) {
      int kPivot=regionIndex[k];
      if(!mark[kPivot]) {
        stack[0]=kPivot;
        CoinBigIndex j=-10;
        next[0]=j;
        int nStack=0;
        while (nStack>=0) {
          /* take off stack */
          if (j>=startR[kPivot]) {
            int jPivot=indexRowR[j--];
            /* put back on stack */
            next[nStack] =j;
            if (!mark[jPivot]) {
            /* and new one */
            kPivot=jPivot;
            j=-10;
            stack[++nStack]=kPivot;
            mark[kPivot]=1;
            next[nStack]=j;
            }
          } else if (j==-10) {
            // before first - see if followon
            int jPivot = permuteBack[kPivot];
            if (jPivot<numberRows_) {
            // no
            j=startR[kPivot]+numberInColumnPlus[kPivot]-1;
            next[nStack]=j;
            } else {
            // add to list
            if (!mark[jPivot]) {
              /* and new one */
              kPivot=jPivot;
              j=-10;
              stack[++nStack]=kPivot;
              mark[kPivot]=1;
              next[nStack]=j;
            } else {
              j=startR[kPivot]+numberInColumnPlus[kPivot]-1;
              next[nStack]=j;
            }
            }
          } else {
            // finished
            list[nList++]=kPivot;
            mark[kPivot]=1;
            --nStack;
            if (nStack>=0) {
            kPivot=stack[nStack];
            j=next[nStack];
            }
          }
        }
      }
      }
      numberNonZero=0;
      for (int i=nList-1;i>=0;i--) {
      int iPivot = list[i];
      mark[iPivot]=0;
      CoinFactorizationDouble pivotValue;
      if (iPivot<numberRows_) {
        pivotValue = region[iPivot];
      } else {
        int before = permute[iPivot];
        pivotValue = region[iPivot] + region[before];
        region[before]=0.0;
      }
      if ( fabs ( pivotValue ) > tolerance ) {
        region[iPivot] = pivotValue;
        CoinBigIndex start = startR[iPivot];
        int number = numberInColumnPlus[iPivot];
        CoinBigIndex end = start + number;
        CoinBigIndex j;
        for (j=start ; j < end; j ++ ) {
          int iRow = indexRowR[j];
          CoinFactorizationDouble value = elementR[j];
          region[iRow] -= value * pivotValue;
        }     
        regionIndex[numberNonZero++] = iPivot;
      } else {
        region[iPivot] = 0.0;
      }       
      }
    }
#else
    {
      
      // use sparse_ as temporary area
      // mark known to be zero
      int * COIN_RESTRICT stack = sparse_.array();  /* pivot */
      int * COIN_RESTRICT list = stack + maximumRowsExtra_;  /* final list */
      CoinBigIndex * COIN_RESTRICT next = reinterpret_cast<CoinBigIndex *> (list + maximumRowsExtra_);  /* jnext */
      char * COIN_RESTRICT mark = reinterpret_cast<char *> (next + maximumRowsExtra_);
      // mark all rows which will be permuted
      for (int i = numberRows_; i < numberRowsExtra_; i++ ) {
      int iRow = permute[i];
      mark[iRow]=1;
      }
      // we have another copy of R in R
      const CoinFactorizationDouble * elementR = elementR_ + lengthAreaR_;
      const int * indexRowR = indexRowR_ + lengthAreaR_;
      const CoinBigIndex * startR = startColumnR_.array()+maximumPivots_+1;
      // For current list order does not matter as
      // only affects end
      int newNumber=0;
      for (int i = 0; i < numberNonZero; i++ ) {
      int iRow = regionIndex[i];
      assert (region[iRow]);
      if (!mark[iRow])
        regionIndex[newNumber++]=iRow;
      int number = numberInColumnPlus[iRow];
      if (number) {
        CoinFactorizationDouble pivotValue = region[iRow];
        CoinBigIndex start=startR[iRow];
        CoinBigIndex end = start+number;
        for (CoinBigIndex j = start; j < end; j ++ ) {
          CoinFactorizationDouble value = elementR[j];
          int jRow = indexRowR[j];
          region[jRow] -= pivotValue*value;
        }
      }
      }
      numberNonZero = newNumber;
      for (int i = numberRows_; i < numberRowsExtra_; i++ ) {
      //move using permute_ (stored in inverse fashion)
      int iRow = permute[i];
      CoinFactorizationDouble pivotValue = region[iRow]+region[i];
      //zero out pre-permuted
      region[iRow] = 0.0;
      if ( fabs ( pivotValue ) > tolerance ) {
        region[i] = pivotValue;
        if (!mark[i])
          regionIndex[numberNonZero++] = i;
        int number = numberInColumnPlus[i];
        CoinBigIndex start=startR[i];
        CoinBigIndex end = start+number;
        for (CoinBigIndex j = start; j < end; j ++ ) {
          CoinFactorizationDouble value = elementR[j];
          int jRow = indexRowR[j];
          region[jRow] -= pivotValue*value;
        }
      } else {
        region[i] = 0.0;
      }
      mark[iRow]=0;
      }
    }
#endif
    break;
  case 1:
    {
      // no sparse region
      // we have another copy of R in R
      const CoinFactorizationDouble * elementR = elementR_ + lengthAreaR_;
      const int * indexRowR = indexRowR_ + lengthAreaR_;
      const CoinBigIndex * startR = startColumnR_.array()+maximumPivots_+1;
      // For current list order does not matter as
      // only affects end
      for (int i = 0; i < numberNonZero; i++ ) {
      int iRow = regionIndex[i];
      assert (region[iRow]);
      int number = numberInColumnPlus[iRow];
      if (number) {
        CoinFactorizationDouble pivotValue = region[iRow];
        CoinBigIndex start=startR[iRow];
        CoinBigIndex end = start+number;
        for (CoinBigIndex j = start; j < end; j ++ ) {
          CoinFactorizationDouble value = elementR[j];
          int jRow = indexRowR[j];
          region[jRow] -= pivotValue*value;
        }
      }
      }
      for (int i = numberRows_; i < numberRowsExtra_; i++ ) {
      //move using permute_ (stored in inverse fashion)
      int iRow = permute[i];
      CoinFactorizationDouble pivotValue = region[iRow]+region[i];
      //zero out pre-permuted
      region[iRow] = 0.0;
      if ( fabs ( pivotValue ) > tolerance ) {
        region[i] = pivotValue;
        regionIndex[numberNonZero++] = i;
        int number = numberInColumnPlus[i];
        CoinBigIndex start=startR[i];
        CoinBigIndex end = start+number;
        for (CoinBigIndex j = start; j < end; j ++ ) {
          CoinFactorizationDouble value = elementR[j];
          int jRow = indexRowR[j];
          region[jRow] -= pivotValue*value;
        }
      } else {
        region[i] = 0.0;
      }
      }
    }
    break;
  case 2:
    {
      CoinBigIndex start = startColumn[numberRows_];
      for (int i = numberRows_; i < numberRowsExtra_; i++ ) {
      //move using permute_ (stored in inverse fashion)
      CoinBigIndex end = startColumn[i+1];
      int iRow = permute[i];
      CoinFactorizationDouble pivotValue = region[iRow];
      //zero out pre-permuted
      region[iRow] = 0.0;

      for (CoinBigIndex j = start; j < end; j ++ ) {
        CoinFactorizationDouble value = element[j];
        int jRow = indexRow[j];
        value *= region[jRow];
        pivotValue -= value;
      }
      start=end;
      if ( fabs ( pivotValue ) > tolerance ) {
        region[i] = pivotValue;
        regionIndex[numberNonZero++] = i;
      } else {
      region[i] = 0.0;
      }
      }
    }
    break;
  }
  if (method) {
    // pack down
    int n=numberNonZero;
    numberNonZero=0;
    for (int i=0;i<n;i++) {
      int indexValue = regionIndex[i];
      double value = region[indexValue];
      if (value) 
      regionIndex[numberNonZero++]=indexValue;
    }
  }
  //set counts
  regionSparse->setNumElements ( numberNonZero );
}
//  updateColumnR.  Updates part of column (FTRANR)
void
CoinFactorization::updateColumnRFT ( CoinIndexedVector * regionSparse,
                             int * COIN_RESTRICT regionIndex) 
{
  double * COIN_RESTRICT region = regionSparse->denseVector (  );
  //int *regionIndex = regionSparse->getIndices (  );
  CoinBigIndex * COIN_RESTRICT startColumnU = startColumnU_.array();
  int numberNonZero = regionSparse->getNumElements (  );

  if ( numberR_ ) {
    double tolerance = zeroTolerance_;
    
    const CoinBigIndex * startColumn = startColumnR_.array()-numberRows_;
    const int * indexRow = indexRowR_;
    const CoinFactorizationDouble * element = elementR_;
    const int * permute = permute_.array();
    

    // Work out very dubious idea of what would be fastest
    int method=-1;
    // Size of R
    double sizeR=startColumnR_.array()[numberR_];
    // Average
    double averageR = sizeR/(static_cast<double> (numberRowsExtra_));
    // weights (relative to actual work)
    double setMark = 0.1; // setting mark
    double test1= 1.0; // starting ftran (without testPivot)
    double testPivot = 2.0; // Seeing if zero etc
    double startDot=2.0; // For starting dot product version
    // For final scan
    double final = numberNonZero*1.0;
    double methodTime[3];
    // For second type
    methodTime[1] = numberPivots_ * (testPivot + ((static_cast<double> (numberNonZero))/(static_cast<double> (numberRows_))
                                      * averageR));
    methodTime[1] += numberNonZero *(test1 + averageR);
    // For first type
    methodTime[0] = methodTime[1] + (numberNonZero+numberPivots_)*setMark;
    methodTime[1] += numberNonZero*final;
    // third
    methodTime[2] = sizeR + numberPivots_*startDot + numberNonZero*final;
    // switch off if necessary
    if (!numberInColumnPlus_.array()) {
      methodTime[0]=1.0e100;
      methodTime[1]=1.0e100;
    } else if (!sparse_.array()) {
      methodTime[0]=1.0e100;
    }
    const int * numberInColumnPlus = numberInColumnPlus_.array();
    int * numberInColumn = numberInColumn_.array();
    // adjust for final scan
    methodTime[1] += final;
    double best=1.0e100;
    for (int i=0;i<3;i++) {
      if (methodTime[i]<best) {
      best=methodTime[i];
      method=i;
      }
    }
    assert (method>=0);
    
    switch (method) {
    case 0:
      {
      // use sparse_ as temporary area
      // mark known to be zero
      int * COIN_RESTRICT stack = sparse_.array();  /* pivot */
      int * COIN_RESTRICT list = stack + maximumRowsExtra_;  /* final list */
      CoinBigIndex * COIN_RESTRICT next = reinterpret_cast<CoinBigIndex *> (list + maximumRowsExtra_);  /* jnext */
      char * COIN_RESTRICT mark = reinterpret_cast<char *> (next + maximumRowsExtra_);
      // mark all rows which will be permuted
      for (int i = numberRows_; i < numberRowsExtra_; i++ ) {
        int iRow = permute[i];
        mark[iRow]=1;
      }
      // we have another copy of R in R
      const CoinFactorizationDouble * elementR = elementR_ + lengthAreaR_;
      const int * indexRowR = indexRowR_ + lengthAreaR_;
      const CoinBigIndex * startR = startColumnR_.array()+maximumPivots_+1;
      //save in U
      //in at end
      int iColumn = numberColumnsExtra_;

      startColumnU[iColumn] = startColumnU[maximumColumnsExtra_];
      CoinBigIndex start = startColumnU[iColumn];
  
      //int * putIndex = indexRowU_ + start;
      CoinFactorizationDouble * COIN_RESTRICT putElement = elementU_.array() + start;
      // For current list order does not matter as
      // only affects end
      int newNumber=0;
      for (int i = 0; i < numberNonZero; i++ ) {
        int iRow = regionIndex[i];
        CoinFactorizationDouble pivotValue = region[iRow];
        assert (region[iRow]);
        if (!mark[iRow]) {
          //putIndex[newNumber]=iRow;
          putElement[newNumber]=pivotValue;;
          regionIndex[newNumber++]=iRow;
        }
        int number = numberInColumnPlus[iRow];
        if (number) {
          CoinBigIndex start=startR[iRow];
          CoinBigIndex end = start+number;
          for (CoinBigIndex j = start; j < end; j ++ ) {
            CoinFactorizationDouble value = elementR[j];
            int jRow = indexRowR[j];
            region[jRow] -= pivotValue*value;
          }
        }
      }
      numberNonZero = newNumber;
      for (int i = numberRows_; i < numberRowsExtra_; i++ ) {
        //move using permute_ (stored in inverse fashion)
        int iRow = permute[i];
        CoinFactorizationDouble pivotValue = region[iRow]+region[i];
        //zero out pre-permuted
        region[iRow] = 0.0;
        if ( fabs ( pivotValue ) > tolerance ) {
          region[i] = pivotValue;
          if (!mark[i]) {
            //putIndex[numberNonZero]=i;
            putElement[numberNonZero]=pivotValue;;
            regionIndex[numberNonZero++]=i;
          }
          int number = numberInColumnPlus[i];
          CoinBigIndex start=startR[i];
          CoinBigIndex end = start+number;
          for (CoinBigIndex j = start; j < end; j ++ ) {
            CoinFactorizationDouble value = elementR[j];
            int jRow = indexRowR[j];
            region[jRow] -= pivotValue*value;
          }
        } else {
          region[i] = 0.0;
        }
        mark[iRow]=0;
      }
      numberInColumn[iColumn] = numberNonZero;
      startColumnU[maximumColumnsExtra_] = start + numberNonZero;
      }
      break;
    case 1:
      {
      // no sparse region
      // we have another copy of R in R
      const CoinFactorizationDouble * elementR = elementR_ + lengthAreaR_;
      const int * indexRowR = indexRowR_ + lengthAreaR_;
      const CoinBigIndex * startR = startColumnR_.array()+maximumPivots_+1;
      // For current list order does not matter as
      // only affects end
      for (int i = 0; i < numberNonZero; i++ ) {
        int iRow = regionIndex[i];
        assert (region[iRow]);
        int number = numberInColumnPlus[iRow];
        if (number) {
          CoinFactorizationDouble pivotValue = region[iRow];
          CoinBigIndex start=startR[iRow];
          CoinBigIndex end = start+number;
          for (CoinBigIndex j = start; j < end; j ++ ) {
            CoinFactorizationDouble value = elementR[j];
            int jRow = indexRowR[j];
            region[jRow] -= pivotValue*value;
          }
        }
      }
      for (int i = numberRows_; i < numberRowsExtra_; i++ ) {
        //move using permute_ (stored in inverse fashion)
        int iRow = permute[i];
        CoinFactorizationDouble pivotValue = region[iRow]+region[i];
        //zero out pre-permuted
        region[iRow] = 0.0;
        if ( fabs ( pivotValue ) > tolerance ) {
          region[i] = pivotValue;
          regionIndex[numberNonZero++] = i;
          int number = numberInColumnPlus[i];
          CoinBigIndex start=startR[i];
          CoinBigIndex end = start+number;
          for (CoinBigIndex j = start; j < end; j ++ ) {
            CoinFactorizationDouble value = elementR[j];
            int jRow = indexRowR[j];
            region[jRow] -= pivotValue*value;
          }
        } else {
          region[i] = 0.0;
        }
      }
      }
      break;
    case 2:
      {
      CoinBigIndex start = startColumn[numberRows_];
      for (int i = numberRows_; i < numberRowsExtra_; i++ ) {
        //move using permute_ (stored in inverse fashion)
        CoinBigIndex end = startColumn[i+1];
        int iRow = permute[i];
        CoinFactorizationDouble pivotValue = region[iRow];
        //zero out pre-permuted
        region[iRow] = 0.0;
        
        for (CoinBigIndex j = start; j < end; j ++ ) {
          CoinFactorizationDouble value = element[j];
          int jRow = indexRow[j];
          value *= region[jRow];
          pivotValue -= value;
        }
        start=end;
        if ( fabs ( pivotValue ) > tolerance ) {
          region[i] = pivotValue;
          regionIndex[numberNonZero++] = i;
        } else {
          region[i] = 0.0;
        }
      }
      }
      break;
    }
    if (method) {
      // pack down
      int n=numberNonZero;
      numberNonZero=0;
      //save in U
      //in at end
      int iColumn = numberColumnsExtra_;
      
      assert(startColumnU[iColumn] == startColumnU[maximumColumnsExtra_]);
      CoinBigIndex start = startColumnU[iColumn];
      
      int * COIN_RESTRICT putIndex = indexRowU_.array() + start;
      CoinFactorizationDouble * COIN_RESTRICT putElement = elementU_.array() + start;
      for (int i=0;i<n;i++) {
      int indexValue = regionIndex[i];
      double value = region[indexValue];
      if (value) {
        putIndex[numberNonZero]=indexValue;
        putElement[numberNonZero]=value;
        regionIndex[numberNonZero++]=indexValue;
      }
      }
      numberInColumn[iColumn] = numberNonZero;
      startColumnU[maximumColumnsExtra_] = start + numberNonZero;
    }
    //set counts
    regionSparse->setNumElements ( numberNonZero );
  } else {
    // No R but we still need to save column
    //save in U
    //in at end
    int * COIN_RESTRICT numberInColumn = numberInColumn_.array();
    numberNonZero = regionSparse->getNumElements (  );
    int iColumn = numberColumnsExtra_;
    
    assert(startColumnU[iColumn] == startColumnU[maximumColumnsExtra_]);
    CoinBigIndex start = startColumnU[iColumn];
    numberInColumn[iColumn] = numberNonZero;
    startColumnU[maximumColumnsExtra_] = start + numberNonZero;
    
    int * COIN_RESTRICT putIndex = indexRowU_.array() + start;
    CoinFactorizationDouble * COIN_RESTRICT putElement = elementU_.array() + start;
    for (int i=0;i<numberNonZero;i++) {
      int indexValue = regionIndex[i];
      double value = region[indexValue];
      putIndex[i]=indexValue;
      putElement[i]=value;
    }
  }
}
/* Updates one column (FTRAN) from region2 and permutes.
   region1 starts as zero
   Note - if regionSparse2 packed on input - will be packed on output
   - returns un-permuted result in region2 and region1 is zero */
02021 int CoinFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
                              CoinIndexedVector * regionSparse2)
{
  //permute and move indices into index array
  int * COIN_RESTRICT regionIndex = regionSparse->getIndices (  );
  int numberNonZero = regionSparse2->getNumElements();
  const int *permute = permute_.array();
  int * COIN_RESTRICT index = regionSparse2->getIndices();
  double * COIN_RESTRICT region = regionSparse->denseVector();
  double * COIN_RESTRICT array = regionSparse2->denseVector();
  CoinBigIndex * COIN_RESTRICT startColumnU = startColumnU_.array();
  bool doFT=doForrestTomlin_;
  // see if room
  if (doFT) {
    int iColumn = numberColumnsExtra_;
    
    startColumnU[iColumn] = startColumnU[maximumColumnsExtra_];
    CoinBigIndex start = startColumnU[iColumn];
    CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
    doFT = space>=0;
    if (doFT) {
      regionIndex = indexRowU_.array() + start;
    } else {
      startColumnU[maximumColumnsExtra_] = lengthAreaU_+1;
    }
  }

#ifndef CLP_FACTORIZATION
  bool packed = regionSparse2->packedMode();
  if (packed) {
#else
    assert (regionSparse2->packedMode());
#endif
    for (int j = 0; j < numberNonZero; j ++ ) {
      int iRow = index[j];
      double value = array[j];
      array[j]=0.0;
      iRow = permute[iRow];
      region[iRow] = value;
      regionIndex[j] = iRow;
    }
#ifndef CLP_FACTORIZATION
  } else {
    for (int j = 0; j < numberNonZero; j ++ ) {
      int iRow = index[j];
      double value = array[iRow];
      array[iRow]=0.0;
      iRow = permute[iRow];
      region[iRow] = value;
      regionIndex[j] = iRow;
    }
  }
#endif
  regionSparse->setNumElements ( numberNonZero );
  if (collectStatistics_) {
    numberFtranCounts_++;
    ftranCountInput_ += numberNonZero;
  }
    
  //  ******* L
#if 0
  {
    double *region = regionSparse->denseVector (  );
    //int *regionIndex = regionSparse->getIndices (  );
    int numberNonZero = regionSparse->getNumElements (  );
    for (int i=0;i<numberNonZero;i++) {
      int iRow = regionIndex[i];
      assert (region[iRow]);
    }
  }
#endif
  updateColumnL ( regionSparse, regionIndex );
#if 0
  {
    double *region = regionSparse->denseVector (  );
    //int *regionIndex = regionSparse->getIndices (  );
    int numberNonZero = regionSparse->getNumElements (  );
    for (int i=0;i<numberNonZero;i++) {
      int iRow = regionIndex[i];
      assert (region[iRow]);
    }
  }
#endif
  if (collectStatistics_) 
    ftranCountAfterL_ += regionSparse->getNumElements();
  //permute extra
  //row bits here
  if ( doFT ) 
    updateColumnRFT ( regionSparse, regionIndex );
  else
    updateColumnR ( regionSparse );
  if (collectStatistics_) 
    ftranCountAfterR_ += regionSparse->getNumElements();
  //  ******* U
  updateColumnU ( regionSparse, regionIndex);
  if (!doForrestTomlin_) {
    // Do PFI after everything else
    updateColumnPFI(regionSparse);
  }
  permuteBack(regionSparse,regionSparse2);
  // will be negative if no room
  if ( doFT ) 
    return regionSparse2->getNumElements();
  else 
    return -regionSparse2->getNumElements();
}

Generated by  Doxygen 1.6.0   Back to index