// 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 part of column (FTRANL) when densish
void
CoinFactorization::updateColumnLDensish ( CoinIndexedVector * regionSparse ,
int * regionIndex)
const
{
double *region = regionSparse->denseVector ( );
int number = regionSparse->getNumElements ( );
int numberNonZero;
double tolerance = zeroTolerance_;
numberNonZero = 0;
int k;
int i , iPivot;
const CoinBigIndex *startColumn = startColumnL_.array();
const int *indexRow = indexRowL_.array();
const double *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 (k=0;k<number;k++) {
iPivot=regionIndex[k];
if (iPivot<baseL_)
regionIndex[numberNonZero++]=iPivot;
else
smallestIndex = CoinMin(iPivot,smallestIndex);
}
// now others
for ( i = smallestIndex; i < last; i++ ) {
double 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];
double result = region[iRow];
double value = element[j];
region[iRow] = result - value * pivotValue;
}
regionIndex[numberNonZero++] = i;
} else {
region[i] = 0.0;
}
}
// and dense
for ( ; i < numberRows_; i++ ) {
double 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 * regionIndex)
const
{
double *region = regionSparse->denseVector ( );
int number = regionSparse->getNumElements ( );
int numberNonZero;
double tolerance = zeroTolerance_;
numberNonZero = 0;
int k;
int i , iPivot;
const CoinBigIndex *startColumn = startColumnL_.array();
const int *indexRow = indexRowL_.array();
const double *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 * mark = (CoinCheckZero *) (sparse_.array()+(2+nInBig)*maximumRowsExtra_);
int smallestIndex = numberRowsExtra_;
// do easy ones
for (k=0;k<number;k++) {
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] |= 1<<iBit;
} else {
mark[iWord] = 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);
for ( i = smallestIndex; i < jLast; i++ ) {
double 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];
double result = region[iRow];
double value = element[j];
region[iRow] = result - value * pivotValue;
int iWord = iRow>>CHECK_SHIFT;
int iBit = iRow-(iWord<<CHECK_SHIFT);
if (mark[iWord]) {
mark[iWord] |= 1<<iBit;
} else {
mark[iWord] = 1<<iBit;
}
}
regionIndex[numberNonZero++] = i;
} else {
region[i] = 0.0;
}
}
int kLast = last>>CHECK_SHIFT;
if (jLast<last) {
// now do in chunks
for (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++ ) {
double 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];
double result = region[iRow];
double value = element[j];
region[iRow] = result - value * pivotValue;
int iWord = iRow>>CHECK_SHIFT;
int iBit = iRow-(iWord<<CHECK_SHIFT);
if (mark[iWord]) {
mark[iWord] |= 1<<iBit;
} else {
mark[iWord] = 1<<iBit;
}
}
regionIndex[numberNonZero++] = i;
} else {
region[i] = 0.0;
}
}
mark[k]=0; // zero out marked
}
}
i = kLast<<CHECK_SHIFT;
}
for ( ; i < last; i++ ) {
double 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];
double result = region[iRow];
double 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 * regionIndex)
const
{
double *region = regionSparse->denseVector ( );
int number = regionSparse->getNumElements ( );
int numberNonZero;
double tolerance = zeroTolerance_;
numberNonZero = 0;
int j, k;
int i;
const CoinBigIndex *startColumn = startColumnL_.array();
const int *indexRow = indexRowL_.array();
const double *element = elementL_.array();
// use sparse_ as temporary area
// mark known to be zero
int * stack = sparse_.array(); /* pivot */
int * list = stack + maximumRowsExtra_; /* final list */
CoinBigIndex * next = (CoinBigIndex *) (list + maximumRowsExtra_); /* jnext */
char * mark = (char *) (next + maximumRowsExtra_);
int nList;
#ifdef COIN_DEBUG
for (i=0;i<maximumRowsExtra_;i++) {
assert (!mark[i]);
}
#endif
int nStack;
nList=0;
for (k=0;k<number;k++) {
int kPivot=regionIndex[k];
if (kPivot>=baseL_) {
if(!mark[kPivot]) {
stack[0]=kPivot;
CoinBigIndex j=startColumn[kPivot+1]-1;
nStack=0;
while (nStack>=0) {
/* take off stack */
if (j>=startColumn[kPivot]) {
int jPivot=indexRow[j--];
assert (jPivot>=baseL_);
/* put back on stack */
next[nStack] =j;
if (!mark[jPivot]) {
/* and new one */
kPivot=jPivot;
j = startColumn[kPivot+1]-1;
stack[++nStack]=kPivot;
mark[kPivot]=1;
next[nStack]=j;
}
} else {
/* finished so mark */
list[nList++]=kPivot;
mark[kPivot]=1;
--nStack;
if (nStack>=0) {
kPivot=stack[nStack];
j=next[nStack];
}
}
}
}
} else {
// just put on list
regionIndex[numberNonZero++]=kPivot;
}
}
for (i=nList-1;i>=0;i--) {
int iPivot = list[i];
mark[iPivot]=0;
double pivotValue = region[iPivot];
if ( fabs ( pivotValue ) > tolerance ) {
regionIndex[numberNonZero++]=iPivot;
for ( j = startColumn[iPivot]; j < startColumn[iPivot+1]; j ++ ) {
int iRow = indexRow[j];
double value = element[j];
region[iRow] -= value * pivotValue;
}
} else {
region[iPivot]=0.0;
}
}
regionSparse->setNumElements ( numberNonZero );
}
// updateColumnL. Updates part of column (FTRANL)
void
CoinFactorization::updateColumnL ( CoinIndexedVector * regionSparse,
int * regionIndex) const
{
if (numberL_) {
int number = regionSparse->getNumElements ( );
int goSparse;
// Guess at number at end
if (sparseThreshold_>0) {
if (ftranAverageAfterL_) {
int newNumber = (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 *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) {
//int iopt=0;
//dges(denseArea_,&numberDense_,&numberDense_,densePermute_,
// ®ion[lastSparse],&iopt);
char trans = 'N';
int ione=1;
int info;
F77_FUNC(dgetrs,DGETRS)(&trans,&numberDense_,&ione,denseArea_,&numberDense_,
densePermute_,region+lastSparse,&numberDense_,&info,1);
for (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
}
int CoinFactorization::checkPivot(double saveFromU,
double oldPivot) const
{
int status;
if ( fabs ( saveFromU ) > 1.0e-8 ) {
double checkTolerance;
if ( numberRowsExtra_ < numberRows_ + 2 ) {
checkTolerance = 1.0e-5;
} else if ( numberRowsExtra_ < numberRows_ + 10 ) {
checkTolerance = 1.0e-6;
} else if ( numberRowsExtra_ < numberRows_ + 50 ) {
checkTolerance = 1.0e-8;
} else {
checkTolerance = 1.0e-10;
}
checkTolerance *= relaxCheck_;
if ( fabs ( 1.0 - fabs ( saveFromU / oldPivot ) ) < checkTolerance ) {
status = 0;
} else {
#if COIN_DEBUG
std::cout <<"inaccurate pivot "<< oldPivot << " "
<< saveFromU << std::endl;
#endif
if ( fabs ( fabs ( oldPivot ) - fabs ( saveFromU ) ) < 1.0e-12 ||
fabs ( 1.0 - fabs ( saveFromU / oldPivot ) ) < 1.0e-8 ) {
status = 1;
} else {
status = 2;
}
}
} else {
//error
status = 2;
#if COIN_DEBUG
std::cout <<"inaccurate pivot "<< saveFromU / oldPivot
<< " " << saveFromU << std::endl;
#endif
}
return status;
}
// replaceColumn. Replaces one Column to basis
// returns 0=OK, 1=Probably OK, 2=singular, 3=no room
//partial update already in U
int
CoinFactorization::replaceColumn ( CoinIndexedVector * regionSparse,
int pivotRow,
double pivotCheck ,
bool checkBeforeModifying)
{
CoinBigIndex * startColumnU = startColumnU_.array();
CoinBigIndex *startColumn;
int *indexRow;
double *element;
//return at once if too many iterations
if ( numberColumnsExtra_ >= maximumColumnsExtra_ ) {
return 5;
}
if ( lengthAreaU_ < startColumnU[maximumColumnsExtra_] ) {
return 3;
}
int * numberInRow = numberInRow_.array();
int * numberInColumn = numberInColumn_.array();
int * numberInColumnPlus = numberInColumnPlus_.array();
int realPivotRow;
realPivotRow = pivotColumn_.array()[pivotRow];
//zeroed out region
double *region = regionSparse->denseVector ( );
element = elementU_.array();
//take out old pivot column
// If we have done no pivots then always check before modification
if (!numberPivots_)
checkBeforeModifying=true;
totalElements_ -= numberInColumn[realPivotRow];
double * pivotRegion = pivotRegion_.array();
double oldPivot = pivotRegion[realPivotRow];
// for accuracy check
pivotCheck = pivotCheck / oldPivot;
#if COIN_DEBUG>1
int checkNumber=1000000;
//if (numberL_) checkNumber=-1;
if (numberR_>=checkNumber) {
printf("pivot row %d, check %g - alpha region:\n",
realPivotRow,pivotCheck);
/*int i;
for (i=0;i<numberRows_;i++) {
if (pivotRegion[i])
printf("%d %g\n",i,pivotRegion[i]);
}*/
}
#endif
pivotRegion[realPivotRow] = 0.0;
CoinBigIndex i;
CoinBigIndex saveEnd = startColumnU[realPivotRow]
+ numberInColumn[realPivotRow];
// not necessary at present - but take no chances for future
numberInColumn[realPivotRow] = 0;
//get entries in row (pivot not stored)
CoinBigIndex *startRow = startRowU_.array();
CoinBigIndex start;
CoinBigIndex end;
start = startRow[realPivotRow];
end = start + numberInRow[realPivotRow];
int numberNonZero = 0;
int *indexColumn = indexColumnU_.array();
CoinBigIndex *convertRowToColumn = convertRowToColumnU_.array();
int *regionIndex = regionSparse->getIndices ( );
#if COIN_DEBUG>1
if (numberR_>=checkNumber)
printf("Before btranu\n");
#endif
int smallestIndex=numberRowsExtra_;
if (!checkBeforeModifying) {
for ( i = start; i < end ; i ++ ) {
int iColumn = indexColumn[i];
smallestIndex = CoinMin(smallestIndex,iColumn);
CoinBigIndex j = convertRowToColumn[i];
region[iColumn] = element[j];
#if COIN_DEBUG>1
if (numberR_>=checkNumber)
printf("%d %g\n",iColumn,region[iColumn]);
#endif
element[j] = 0.0;
regionIndex[numberNonZero++] = iColumn;
}
} else {
for ( i = start; i < end ; i ++ ) {
int iColumn = indexColumn[i];
smallestIndex = CoinMin(smallestIndex,iColumn);
CoinBigIndex j = convertRowToColumn[i];
region[iColumn] = element[j];
#if COIN_DEBUG>1
if (numberR_>=checkNumber)
printf("%d %g\n",iColumn,region[iColumn]);
#endif
regionIndex[numberNonZero++] = iColumn;
}
}
//do BTRAN - finding first one to use
regionSparse->setNumElements ( numberNonZero );
updateColumnTransposeU ( regionSparse, smallestIndex );
numberNonZero = regionSparse->getNumElements ( );
double saveFromU = 0.0;
CoinBigIndex startU = startColumnU[numberColumnsExtra_];
int *indexU = &indexRowU_.array()[startU];
double *elementU = &elementU_.array()[startU];
// Do accuracy test here if caller is paranoid
if (checkBeforeModifying) {
double tolerance = zeroTolerance_;
int number = numberInColumn[numberColumnsExtra_];
for ( i = 0; i < number; i++ ) {
int iRow = indexU[i];
//if (numberCompressions_==99&&lengthU_==278)
//printf("row %d saveFromU %g element %g region %g\n",
// iRow,saveFromU,elementU[i],region[iRow]);
if ( fabs ( elementU[i] ) > tolerance ) {
if ( iRow != realPivotRow ) {
saveFromU -= elementU[i] * region[iRow];
} else {
saveFromU += elementU[i];
}
}
}
//check accuracy
int status = checkPivot(saveFromU,pivotCheck);
if (status) {
// restore some things
pivotRegion[realPivotRow] = oldPivot;
number = saveEnd-startColumnU[realPivotRow];
totalElements_ += number;
numberInColumn[realPivotRow]=number;
regionSparse->clear();
return status;
} else {
// do what we would have done by now
for ( i = start; i < end ; i ++ ) {
CoinBigIndex j = convertRowToColumn[i];
element[j] = 0.0;
}
}
}
// Now zero out column of U
//take out old pivot column
for ( i = startColumnU[realPivotRow]; i < saveEnd ; i ++ ) {
element[i] = 0.0;
}
//zero out pivot Row (before or after?)
//add to R
startColumn = startColumnR_.array();
indexRow = indexRowR_;
element = elementR_;
CoinBigIndex l = lengthR_;
int number = numberR_;
startColumn[number] = l; //for luck and first time
number++;
startColumn[number] = l + numberNonZero;
numberR_ = number;
lengthR_ = l + numberNonZero;
totalElements_ += numberNonZero;
if ( lengthR_ >= lengthAreaR_ ) {
//not enough room
regionSparse->clear();
return 3;
}
#if COIN_DEBUG>1
if (numberR_>=checkNumber)
printf("After btranu\n");
#endif
for ( i = 0; i < numberNonZero; i++ ) {
int iRow = regionIndex[i];
#if COIN_DEBUG>1
if (numberR_>=checkNumber)
printf("%d %g\n",iRow,region[iRow]);
#endif
indexRow[l] = iRow;
element[l] = region[iRow];
l++;
}
//take out row
int * nextRow = nextRow_.array();
int * lastRow = lastRow_.array();
int next = nextRow[realPivotRow];
int last = lastRow[realPivotRow];
nextRow[last] = next;
lastRow[next] = last;
numberInRow[realPivotRow]=0;
#if COIN_DEBUG
nextRow[realPivotRow] = 777777;
lastRow[realPivotRow] = 777777;
#endif
//do permute
permute_.array()[numberRowsExtra_] = realPivotRow;
// and other way
permuteBack_.array()[realPivotRow] = numberRowsExtra_;
permuteBack_.array()[numberRowsExtra_] = -1;;
//and for safety
permute_.array()[numberRowsExtra_ + 1] = 0;
pivotColumn_.array()[pivotRow] = numberRowsExtra_;
pivotColumnBack_.array()[numberRowsExtra_] = pivotRow;
startColumn = startColumnU;
indexRow = indexRowU_.array();
element = elementU_.array();
numberU_++;
number = numberInColumn[numberColumnsExtra_];
totalElements_ += number;
lengthU_ += number;
if ( lengthU_ >= lengthAreaU_ ) {
//not enough room
regionSparse->clear();
return 3;
}
saveFromU = 0.0;
//put in pivot
//add row counts
double tolerance = zeroTolerance_;
#if COIN_DEBUG>1
if (numberR_>=checkNumber)
printf("On U\n");
#endif
for ( i = 0; i < number; i++ ) {
int iRow = indexU[i];
#if COIN_DEBUG>1
if (numberR_>=checkNumber)
printf("%d %g\n",iRow,elementU[i]);
#endif
if ( fabs ( elementU[i] ) > tolerance ) {
if ( iRow != realPivotRow ) {
int next = nextRow[iRow];
int iNumberInRow = numberInRow[iRow];
CoinBigIndex space;
CoinBigIndex put = startRow[iRow] + iNumberInRow;
space = startRow[next] - put;
if ( space <= 0 ) {
getRowSpaceIterate ( iRow, iNumberInRow + 4 );
put = startRow[iRow] + iNumberInRow;
}
indexColumn[put] = numberColumnsExtra_;
convertRowToColumn[put] = i + startU;
numberInRow[iRow] = iNumberInRow + 1;
saveFromU = saveFromU - elementU[i] * region[iRow];
} else {
//zero out and save
saveFromU += elementU[i];
elementU[i] = 0.0;
}
} else {
elementU[i] = 0.0;
}
}
//in at end
last = lastRow[maximumRowsExtra_];
nextRow[last] = numberRowsExtra_;
lastRow[maximumRowsExtra_] = numberRowsExtra_;
lastRow[numberRowsExtra_] = last;
nextRow[numberRowsExtra_] = maximumRowsExtra_;
startRow[numberRowsExtra_] = startRow[maximumRowsExtra_];
numberInRow[numberRowsExtra_] = 0;
//column in at beginning (as empty)
int * nextColumn = nextColumn_.array();
int * lastColumn = lastColumn_.array();
next = nextColumn[maximumColumnsExtra_];
lastColumn[next] = numberColumnsExtra_;
nextColumn[maximumColumnsExtra_] = numberColumnsExtra_;
nextColumn[numberColumnsExtra_] = next;
lastColumn[numberColumnsExtra_] = maximumColumnsExtra_;
//check accuracy - but not if already checked (optimization problem)
int status = (checkBeforeModifying) ? 0 : checkPivot(saveFromU,pivotCheck);
if (status!=2) {
double pivotValue = 1.0 / saveFromU;
pivotRegion[numberRowsExtra_] = pivotValue;
//modify by pivot
for ( i = 0; i < number; i++ ) {
elementU[i] *= pivotValue;
}
maximumU_ = CoinMax(maximumU_,startU+number);
numberRowsExtra_++;
numberColumnsExtra_++;
numberGoodU_++;
numberPivots_++;
}
if ( numberRowsExtra_ > numberRows_ + 50 ) {
CoinBigIndex extra = factorElements_ >> 1;
if ( numberRowsExtra_ > numberRows_ + 100 + numberRows_ / 500 ) {
if ( extra < 2 * numberRows_ ) {
extra = 2 * numberRows_;
}
} else {
if ( extra < 5 * numberRows_ ) {
extra = 5 * numberRows_;
}
}
CoinBigIndex added = totalElements_ - factorElements_;
if ( added > extra && added > ( factorElements_ ) << 1 && !status
&& 3*totalElements_ > 2*(lengthAreaU_+lengthAreaL_)) {
status = 3;
if ( messageLevel_ & 4 ) {
std::cout << "Factorization has "<< totalElements_
<< ", basis had " << factorElements_ <<std::endl;
}
}
}
if (numberInColumnPlus&&status<2) {
// we are going to put another copy of R in R
double * elementR = elementR_ + lengthAreaR_;
int * indexRowR = indexRowR_ + lengthAreaR_;
CoinBigIndex * startR = startColumnR_.array()+maximumPivots_+1;
int pivotRow = numberRowsExtra_-1;
for ( i = 0; i < numberNonZero; i++ ) {
int iRow = regionIndex[i];
assert (pivotRow>iRow);
next = nextColumn[iRow];
CoinBigIndex space;
if (next!=maximumColumnsExtra_)
space = startR[next]-startR[iRow];
else
space = lengthAreaR_-startR[iRow];
int numberInR = numberInColumnPlus[iRow];
if (space>numberInR) {
// there is space
CoinBigIndex put=startR[iRow]+numberInR;
numberInColumnPlus[iRow]=numberInR+1;
indexRowR[put]=pivotRow;
elementR[put]=region[iRow];
//add 4 for luck
if (next==maximumColumnsExtra_)
startR[maximumColumnsExtra_] = CoinMin((CoinBigIndex) (put + 4) ,lengthAreaR_);
} else {
// no space - do we shuffle?
if (!getColumnSpaceIterateR(iRow,region[iRow],pivotRow)) {
// printf("Need more space for R\n");
numberInColumnPlus_.conditionalDelete();
regionSparse->clear();
break;
}
}
region[iRow]=0.0;
}
regionSparse->setNumElements(0);
} else {
regionSparse->clear();
}
return status;
}
// updateColumnTranspose. Updates one column transpose (BTRAN)
int
CoinFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
CoinIndexedVector * regionSparse2 )
const
{
//zero region
regionSparse->clear ( );
double *region = regionSparse->denseVector ( );
double * vector = regionSparse2->denseVector();
int * index = regionSparse2->getIndices();
int numberNonZero = regionSparse2->getNumElements();
int i;
const int * pivotColumn = pivotColumn_.array();
//move indices into index array
int *regionIndex = regionSparse->getIndices ( );
int iRow;
bool packed = regionSparse2->packedMode();
if (packed) {
for ( i = 0; i < numberNonZero; i ++ ) {
iRow = index[i];
double value = vector[i];
iRow=pivotColumn[iRow];
vector[i]=0.0;
region[iRow] = value;
regionIndex[i] = iRow;
}
} else {
for ( i = 0; i < numberNonZero; i ++ ) {
iRow = index[i];
double value = vector[iRow];
vector[iRow]=0.0;
iRow=pivotColumn[iRow];
region[iRow] = value;
regionIndex[i] = iRow;
}
}
regionSparse->setNumElements ( numberNonZero );
if (collectStatistics_) {
numberBtranCounts_++;
btranCountInput_ += (double) numberNonZero;
}
if (!doForrestTomlin_) {
// Do PFI before everything else
updateColumnTransposePFI(regionSparse);
numberNonZero = regionSparse->getNumElements();
}
// ******* U
// Apply pivot region - could be combined for speed
int j;
double *pivotRegion = pivotRegion_.array();
int smallestIndex=numberRowsExtra_;
for ( j = 0; j < numberNonZero; j++ ) {
int iRow = regionIndex[j];
smallestIndex = CoinMin(smallestIndex,iRow);
region[iRow] *= pivotRegion[iRow];
}
updateColumnTransposeU ( regionSparse,smallestIndex );
if (collectStatistics_)
btranCountAfterU_ += (double) regionSparse->getNumElements();
//permute extra
//row bits here
updateColumnTransposeR ( regionSparse );
// ******* L
updateColumnTransposeL ( regionSparse );
numberNonZero = regionSparse->getNumElements ( );
if (collectStatistics_)
btranCountAfterL_ += (double) numberNonZero;
const int * permuteBack = pivotColumnBack_.array();
int number=0;
if (packed) {
for (i=0;i<numberNonZero;i++) {
int iRow=regionIndex[i];
double value = region[iRow];
region[iRow]=0.0;
if (fabs(value)>zeroTolerance_) {
iRow=permuteBack[iRow];
vector[number]=value;
index[number++]=iRow;
}
}
} else {
for (i=0;i<numberNonZero;i++) {
int iRow=regionIndex[i];
double value = region[iRow];
region[iRow]=0.0;
if (fabs(value)>zeroTolerance_) {
iRow=permuteBack[iRow];
vector[iRow]=value;
index[number++]=iRow;
}
}
}
regionSparse->setNumElements(0);
regionSparse2->setNumElements(number);
#ifdef COIN_DEBUG
for (i=0;i<numberRowsExtra_;i++) {
assert (!region[i]);
}
#endif
return number;
}
/* Updates part of column transpose (BTRANU) when densish,
assumes index is sorted i.e. region is correct */
void
CoinFactorization::updateColumnTransposeUDensish
( CoinIndexedVector * regionSparse,
int smallestIndex) const
{
double *region = regionSparse->denseVector ( );
int numberNonZero = regionSparse->getNumElements ( );
double tolerance = zeroTolerance_;
int *regionIndex = regionSparse->getIndices ( );
int i,j;
const CoinBigIndex *startRow = startRowU_.array();
const CoinBigIndex *convertRowToColumn = convertRowToColumnU_.array();
const int *indexColumn = indexColumnU_.array();
const double * element = elementU_.array();
int last = numberU_;
double pivotValue;
const int *numberInRow = numberInRow_.array();
numberNonZero = 0;
for (i=smallestIndex ; i < last; i++ ) {
pivotValue = region[i];
if ( fabs ( pivotValue ) > tolerance ) {
CoinBigIndex start = startRow[i];
int numberIn = numberInRow[i];
CoinBigIndex end = start + numberIn;
for (j = start ; j < end; j ++ ) {
int iRow = indexColumn[j];
CoinBigIndex getElement = convertRowToColumn[j];
double value = element[getElement];
region[iRow] -= value * pivotValue;
}
regionIndex[numberNonZero++] = i;
} else {
region[i] = 0.0;
}
}
//set counts
regionSparse->setNumElements ( numberNonZero );
}
/* Updates part of column transpose (BTRANU) when sparsish,
assumes index is sorted i.e. region is correct */
void
CoinFactorization::updateColumnTransposeUSparsish
( CoinIndexedVector * regionSparse,
int smallestIndex) const
{
double *region = regionSparse->denseVector ( );
int numberNonZero = regionSparse->getNumElements ( );
double tolerance = zeroTolerance_;
int *regionIndex = regionSparse->getIndices ( );
int i,j;
const CoinBigIndex *startRow = startRowU_.array();
const CoinBigIndex *convertRowToColumn = convertRowToColumnU_.array();
const int *indexColumn = indexColumnU_.array();
const double * element = elementU_.array();
int last = numberU_;
double pivotValue;
const int *numberInRow = numberInRow_.array();
// mark known to be zero
int nInBig = sizeof(CoinBigIndex)/sizeof(int);
CoinCheckZero * mark = (CoinCheckZero *) (sparse_.array()+(2+nInBig)*maximumRowsExtra_);
for (i=0;i<numberNonZero;i++) {
int iPivot=regionIndex[i];
int iWord = iPivot>>CHECK_SHIFT;
int iBit = iPivot-(iWord<<CHECK_SHIFT);
if (mark[iWord]) {
mark[iWord] |= 1<<iBit;
} else {
mark[iWord] = 1<<iBit;
}
}
numberNonZero = 0;
// Find convenient power of 2
smallestIndex = smallestIndex >> CHECK_SHIFT;
int kLast = last>>CHECK_SHIFT;
// do in chunks
int k;
for (k=smallestIndex;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++ ) {
pivotValue = region[i];
if ( fabs ( pivotValue ) > tolerance ) {
CoinBigIndex start = startRow[i];
int numberIn = numberInRow[i];
CoinBigIndex end = start + numberIn;
for (j = start ; j < end; j ++ ) {
int iRow = indexColumn[j];
CoinBigIndex getElement = convertRowToColumn[j];
double value = element[getElement];
int iWord = iRow>>CHECK_SHIFT;
int iBit = iRow-(iWord<<CHECK_SHIFT);
if (mark[iWord]) {
mark[iWord] |= 1<<iBit;
} else {
mark[iWord] = 1<<iBit;
}
region[iRow] -= value * pivotValue;
}
regionIndex[numberNonZero++] = i;
} else {
region[i] = 0.0;
}
}
mark[k]=0;
}
}
i = kLast<<CHECK_SHIFT;
mark[kLast]=0;
for (; i < last; i++ ) {
pivotValue = region[i];
if ( fabs ( pivotValue ) > tolerance ) {
CoinBigIndex start = startRow[i];
int numberIn = numberInRow[i];
CoinBigIndex end = start + numberIn;
for (j = start ; j < end; j ++ ) {
int iRow = indexColumn[j];
CoinBigIndex getElement = convertRowToColumn[j];
double value = element[getElement];
region[iRow] -= value * pivotValue;
}
regionIndex[numberNonZero++] = i;
} else {
region[i] = 0.0;
}
}
#ifdef COIN_DEBUG
for (i=0;i<maximumRowsExtra_;i++) {
assert (!mark[i]);
}
#endif
//set counts
regionSparse->setNumElements ( numberNonZero );
}
/* Updates part of column transpose (BTRANU) when sparse,
assumes index is sorted i.e. region is correct */
void
CoinFactorization::updateColumnTransposeUSparse (
CoinIndexedVector * regionSparse) const
{
double *region = regionSparse->denseVector ( );
int numberNonZero = regionSparse->getNumElements ( );
double tolerance = zeroTolerance_;
int *regionIndex = regionSparse->getIndices ( );
int i;
const CoinBigIndex *startRow = startRowU_.array();
const CoinBigIndex *convertRowToColumn = convertRowToColumnU_.array();
const int *indexColumn = indexColumnU_.array();
const double * element = elementU_.array();
double pivotValue;
int *numberInRow = numberInRow_.array();
// use sparse_ as temporary area
// mark known to be zero
int * stack = sparse_.array(); /* pivot */
int * list = stack + maximumRowsExtra_; /* final list */
CoinBigIndex * next = (CoinBigIndex *) (list + maximumRowsExtra_); /* jnext */
char * mark = (char *) (next + maximumRowsExtra_);
int nList;
int iPivot;
#ifdef COIN_DEBUG
for (i=0;i<maximumRowsExtra_;i++) {
assert (!mark[i]);
}
#endif
#if 0
{
int i;
for (i=0;i<numberRowsExtra_;i++) {
CoinBigIndex krs = startRow[i];
CoinBigIndex kre = krs + numberInRow[i];
CoinBigIndex k;
for (k=krs;k<kre;k++)
assert (indexColumn[k]>i);
}
}
#endif
int k,nStack;
nList=0;
for (k=0;k<numberNonZero;k++) {
int kPivot=regionIndex[k];
stack[0]=kPivot;
CoinBigIndex j=startRow[kPivot]+numberInRow[kPivot]-1;
next[0]=j;
nStack=1;
while (nStack) {
/* take off stack */
kPivot=stack[--nStack];
if (mark[kPivot]!=1) {
j=next[nStack];
if (j>=startRow[kPivot]) {
kPivot=indexColumn[j--];
/* put back on stack */
next[nStack++] =j;
if (!mark[kPivot]) {
/* and new one */
j=startRow[kPivot]+numberInRow[kPivot]-1;
stack[nStack]=kPivot;
mark[kPivot]=2;
next[nStack++]=j;
}
} else {
// finished
list[nList++]=kPivot;
mark[kPivot]=1;
}
}
}
}
numberNonZero=0;
for (i=nList-1;i>=0;i--) {
iPivot = list[i];
mark[iPivot]=0;
pivotValue = region[iPivot];
if ( fabs ( pivotValue ) > tolerance ) {
CoinBigIndex start = startRow[iPivot];
int numberIn = numberInRow[iPivot];
CoinBigIndex end = start + numberIn;
CoinBigIndex j;
for (j=start ; j < end; j ++ ) {
int iRow = indexColumn[j];
CoinBigIndex getElement = convertRowToColumn[j];
double value = element[getElement];
region[iRow] -= value * pivotValue;
}
regionIndex[numberNonZero++] = iPivot;
} else {
region[iPivot] = 0.0;
}
}
//set counts
regionSparse->setNumElements ( numberNonZero );
}
// updateColumnTransposeU. Updates part of column transpose (BTRANU)
//assumes index is sorted i.e. region is correct
//does not sort by sign
void
CoinFactorization::updateColumnTransposeU ( CoinIndexedVector * regionSparse,
int smallestIndex) const
{
int number = regionSparse->getNumElements ( );
int goSparse;
// Guess at number at end
if (sparseThreshold_>0) {
if (btranAverageAfterU_) {
int newNumber = (int) (number*btranAverageAfterU_);
if (newNumber< sparseThreshold_)
goSparse = 2;
else if (newNumber< sparseThreshold2_)
goSparse = 1;
else
goSparse = 0;
} else {
if (number<sparseThreshold_)
goSparse = 2;
else
goSparse = 0;
}
} else {
goSparse=0;
}
switch (goSparse) {
case 0: // densish
updateColumnTransposeUDensish(regionSparse,smallestIndex);
break;
case 1: // middling
updateColumnTransposeUSparsish(regionSparse,smallestIndex);
break;
case 2: // sparse
updateColumnTransposeUSparse(regionSparse);
break;
}
}
/* updateColumnTransposeLDensish.
Updates part of column transpose (BTRANL) dense by column */
void
CoinFactorization::updateColumnTransposeLDensish
( CoinIndexedVector * regionSparse ) const
{
double *region = regionSparse->denseVector ( );
int *regionIndex = regionSparse->getIndices ( );
int numberNonZero;
double tolerance = zeroTolerance_;
int base;
int first = -1;
numberNonZero=0;
//scan
for (first=numberRows_-1;first>=0;first--) {
if (region[first])
break;
}
if ( first >= 0 ) {
base = baseL_;
const CoinBigIndex *startColumn = startColumnL_.array();
const int *indexRow = indexRowL_.array();
const double *element = elementL_.array();
int last = baseL_ + numberL_;
if ( first >= last ) {
first = last - 1;
}
int i;
double pivotValue;
for (i = first ; i >= base; i-- ) {
CoinBigIndex j;
pivotValue = region[i];
for ( j= startColumn[i] ; j < startColumn[i+1]; j++ ) {
int iRow = indexRow[j];
double value = element[j];
pivotValue -= value * region[iRow];
}
if ( fabs ( pivotValue ) > tolerance ) {
region[i] = pivotValue;
regionIndex[numberNonZero++] = i;
} else {
region[i] = 0.0;
}
}
//may have stopped early
if ( first < base ) {
base = first + 1;
}
if (base > 5) {
i=base-1;
pivotValue=region[i];
bool store = fabs(pivotValue) > tolerance;
for (; i > 0; i-- ) {
bool oldStore = store;
double oldValue = pivotValue;
pivotValue = region[i-1];
store = fabs ( pivotValue ) > tolerance ;
if (oldStore) {
region[i] = oldValue;
regionIndex[numberNonZero++] = i;
} else {
region[i] = 0.0;
}
}
if (store) {
region[0] = pivotValue;
regionIndex[numberNonZero++] = 0;
} else {
region[0] = 0.0;
}
} else {
for (i = base -1 ; i >= 0; i-- ) {
pivotValue = region[i];
if ( fabs ( pivotValue ) > tolerance ) {
region[i] = pivotValue;
regionIndex[numberNonZero++] = i;
} else {
region[i] = 0.0;
}
}
}
}
//set counts
regionSparse->setNumElements ( numberNonZero );
}
/* updateColumnTransposeLByRow.
Updates part of column transpose (BTRANL) densish but by row */
void
CoinFactorization::updateColumnTransposeLByRow
( CoinIndexedVector * regionSparse ) const
{
double *region = regionSparse->denseVector ( );
int *regionIndex = regionSparse->getIndices ( );
int numberNonZero;
double tolerance = zeroTolerance_;
int first = -1;
// use row copy of L
const double * element = elementByRowL_.array();
const CoinBigIndex * startRow = startRowL_.array();
const int * column = indexColumnL_.array();
int i;
CoinBigIndex j;
for (first=numberRows_-1;first>=0;first--) {
if (region[first])
break;
}
numberNonZero=0;
for (i=first;i>=0;i--) {
double pivotValue = region[i];
if ( fabs ( pivotValue ) > tolerance ) {
regionIndex[numberNonZero++] = i;
for (j = startRow[i + 1]-1;j >= startRow[i]; j--) {
int iRow = column[j];
double value = element[j];
region[iRow] -= pivotValue*value;
}
} else {
region[i] = 0.0;
}
}
//set counts
regionSparse->setNumElements ( numberNonZero );
}
// Updates part of column transpose (BTRANL) when sparsish by row
void
CoinFactorization::updateColumnTransposeLSparsish
( CoinIndexedVector * regionSparse ) const
{
double *region = regionSparse->denseVector ( );
int *regionIndex = regionSparse->getIndices ( );
int numberNonZero = regionSparse->getNumElements();
double tolerance = zeroTolerance_;
// use row copy of L
const double * element = elementByRowL_.array();
const CoinBigIndex * startRow = startRowL_.array();
const int * column = indexColumnL_.array();
int i;
CoinBigIndex j;
// mark known to be zero
int nInBig = sizeof(CoinBigIndex)/sizeof(int);
CoinCheckZero * mark = (CoinCheckZero *) (sparse_.array()+(2+nInBig)*maximumRowsExtra_);
for (i=0;i<numberNonZero;i++) {
int iPivot=regionIndex[i];
int iWord = iPivot>>CHECK_SHIFT;
int iBit = iPivot-(iWord<<CHECK_SHIFT);
if (mark[iWord]) {
mark[iWord] |= 1<<iBit;
} else {
mark[iWord] = 1<<iBit;
}
}
numberNonZero = 0;
// First do down to convenient power of 2
int jLast = (numberRows_-1)>>CHECK_SHIFT;
jLast = (jLast<<CHECK_SHIFT);
for (i=numberRows_-1;i>=jLast;i--) {
double pivotValue = region[i];
if ( fabs ( pivotValue ) > tolerance ) {
regionIndex[numberNonZero++] = i;
for (j = startRow[i + 1]-1;j >= startRow[i]; j--) {
int iRow = column[j];
double value = element[j];
int iWord = iRow>>CHECK_SHIFT;
int iBit = iRow-(iWord<<CHECK_SHIFT);
if (mark[iWord]) {
mark[iWord] |= 1<<iBit;
} else {
mark[iWord] = 1<<iBit;
}
region[iRow] -= pivotValue*value;
}
} else {
region[i] = 0.0;
}
}
// and in chunks
jLast = jLast>>CHECK_SHIFT;
int k ;
for (k=jLast-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 pivotValue = region[i];
if ( fabs ( pivotValue ) > tolerance ) {
regionIndex[numberNonZero++] = i;
for (j = startRow[i + 1]-1;j >= startRow[i]; j--) {
int iRow = column[j];
double value = element[j];
int iWord = iRow>>CHECK_SHIFT;
int iBit = iRow-(iWord<<CHECK_SHIFT);
if (mark[iWord]) {
mark[iWord] |= 1<<iBit;
} else {
mark[iWord] = 1<<iBit;
}
region[iRow] -= pivotValue*value;
}
} else {
region[i] = 0.0;
}
}
mark[k]=0;
}
}
mark[jLast]=0;
#ifdef COIN_DEBUG
for (i=0;i<maximumRowsExtra_;i++) {
assert (!mark[i]);
}
#endif
//set counts
regionSparse->setNumElements ( numberNonZero );
}
/* updateColumnTransposeLSparse.
Updates part of column transpose (BTRANL) sparse */
void
CoinFactorization::updateColumnTransposeLSparse
( CoinIndexedVector * regionSparse ) const
{
double *region = regionSparse->denseVector ( );
int *regionIndex = regionSparse->getIndices ( );
int numberNonZero = regionSparse->getNumElements ( );
double tolerance = zeroTolerance_;
// use row copy of L
const double * element = elementByRowL_.array();
const CoinBigIndex * startRow = startRowL_.array();
const int * column = indexColumnL_.array();
int i;
CoinBigIndex j;
// use sparse_ as temporary area
// mark known to be zero
int * stack = sparse_.array(); /* pivot */
int * list = stack + maximumRowsExtra_; /* final list */
CoinBigIndex * next = (CoinBigIndex *) (list + maximumRowsExtra_); /* jnext */
char * mark = (char *) (next + maximumRowsExtra_);
int nList;
int number = numberNonZero;
int k, iPivot;
#ifdef COIN_DEBUG
for (i=0;i<maximumRowsExtra_;i++) {
assert (!mark[i]);
}
#endif
int nStack;
nList=0;
for (k=0;k<number;k++) {
int kPivot=regionIndex[k];
if(!mark[kPivot]&®ion[kPivot]) {
stack[0]=kPivot;
CoinBigIndex j=startRow[kPivot+1]-1;
nStack=0;
while (nStack>=0) {
/* take off stack */
if (j>=startRow[kPivot]) {
int jPivot=column[j--];
/* put back on stack */
next[nStack] =j;
if (!mark[jPivot]) {
/* and new one */
kPivot=jPivot;
j = startRow[kPivot+1]-1;
stack[++nStack]=kPivot;
mark[kPivot]=1;
next[nStack]=j;
}
} else {
/* finished so mark */
list[nList++]=kPivot;
mark[kPivot]=1;
--nStack;
if (nStack>=0) {
kPivot=stack[nStack];
j=next[nStack];
}
}
}
}
}
numberNonZero=0;
for (i=nList-1;i>=0;i--) {
iPivot = list[i];
mark[iPivot]=0;
double pivotValue = region[iPivot];
if ( fabs ( pivotValue ) > tolerance ) {
regionIndex[numberNonZero++] = iPivot;
for ( j = startRow[iPivot]; j < startRow[iPivot+1]; j ++ ) {
int iRow = column[j];
double value = element[j];
region[iRow] -= value * pivotValue;
}
} else {
region[iPivot]=0.0;
}
}
//set counts
regionSparse->setNumElements ( numberNonZero );
}
// updateColumnTransposeL. Updates part of column transpose (BTRANL)
void
CoinFactorization::updateColumnTransposeL ( CoinIndexedVector * regionSparse ) const
{
int number = regionSparse->getNumElements ( );
if (!numberL_&&!numberDense_) {
if (sparse_.array()||number<numberRows_)
return;
}
int goSparse;
// Guess at number at end
// we may need to rethink on dense
if (sparseThreshold_>0) {
if (btranAverageAfterL_) {
int newNumber = (int) (number*btranAverageAfterL_);
if (newNumber< sparseThreshold_)
goSparse = 2;
else if (newNumber< sparseThreshold2_)
goSparse = 1;
else
goSparse = 0;
} else {
if (number<sparseThreshold_)
goSparse = 2;
else
goSparse = 0;
}
} else {
goSparse=-1;
}
#ifdef DENSE_CODE
if (numberDense_) {
//take off list
int lastSparse = numberRows_-numberDense_;
int number = regionSparse->getNumElements();
double *region = regionSparse->denseVector ( );
int *regionIndex = regionSparse->getIndices ( );
int i=0;
bool doDense=false;
if (number<=numberRows_) {
while (i<number) {
int iRow = regionIndex[i];
if (iRow>=lastSparse) {
doDense=true;
regionIndex[i] = regionIndex[--number];
} else {
i++;
}
}
} else {
for (i=numberRows_-1;i>=lastSparse;i--) {
if (region[i]) {
doDense=true;
// numbers are all wrong - do a scan
regionSparse->setNumElements(0);
regionSparse->scan(0,lastSparse,zeroTolerance_);
number=regionSparse->getNumElements();
break;
}
}
if (sparseThreshold_)
goSparse=0;
else
goSparse=-1;
}
if (doDense) {
regionSparse->setNumElements(number);
char trans = 'T';
int ione=1;
int info;
F77_FUNC(dgetrs,DGETRS)(&trans,&numberDense_,&ione,denseArea_,&numberDense_,
densePermute_,region+lastSparse,&numberDense_,&info,1);
//and scan again
if (goSparse>0||!numberL_)
regionSparse->scan(lastSparse,numberRows_,zeroTolerance_);
}
if (!numberL_)
return;
}
#endif
switch (goSparse) {
case -1: // No row copy
updateColumnTransposeLDensish(regionSparse);
break;
case 0: // densish but by row
updateColumnTransposeLByRow(regionSparse);
break;
case 1: // middling(and by row)
updateColumnTransposeLSparsish(regionSparse);
break;
case 2: // sparse
updateColumnTransposeLSparse(regionSparse);
break;
}
}
// getRowSpaceIterate. Gets space for one Row with given length
//may have to do compression (returns true)
//also moves existing vector
bool
CoinFactorization::getRowSpaceIterate ( int iRow,
int extraNeeded )
{
const int * numberInRow = numberInRow_.array();
int number = numberInRow[iRow];
CoinBigIndex *startRow = startRowU_.array();
int *indexColumn = indexColumnU_.array();
CoinBigIndex *convertRowToColumn = convertRowToColumnU_.array();
CoinBigIndex space = lengthAreaU_ - startRow[maximumRowsExtra_];
int * nextRow = nextRow_.array();
int * lastRow = lastRow_.array();
if ( space < extraNeeded + number + 2 ) {
//compression
int iRow = nextRow[maximumRowsExtra_];
CoinBigIndex put = 0;
while ( iRow != maximumRowsExtra_ ) {
//move
CoinBigIndex get = startRow[iRow];
CoinBigIndex getEnd = startRow[iRow] + numberInRow[iRow];
startRow[iRow] = put;
CoinBigIndex i;
for ( i = get; i < getEnd; i++ ) {
indexColumn[put] = indexColumn[i];
convertRowToColumn[put] = convertRowToColumn[i];
put++;
}
iRow = nextRow[iRow];
} /* endwhile */
numberCompressions_++;
startRow[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 = startRow[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 = startRow[iRow];
int * indexColumnU = indexColumnU_.array();
startRow[iRow] = put;
while ( number ) {
number--;
indexColumnU[put] = indexColumnU[get];
convertRowToColumn[put] = convertRowToColumn[get];
put++;
get++;
} /* endwhile */
//add four for luck
startRow[maximumRowsExtra_] = put + extraNeeded + 4;
return true;
}
// getColumnSpaceIterateR. Gets space for one extra R element in Column
//may have to do compression (returns true)
//also moves existing vector
bool
CoinFactorization::getColumnSpaceIterateR ( int iColumn, double value,
int iRow)
{
double * elementR = elementR_ + lengthAreaR_;
int * indexRowR = indexRowR_ + lengthAreaR_;
CoinBigIndex * startR = startColumnR_.array()+maximumPivots_+1;
int * numberInColumnPlus = numberInColumnPlus_.array();
int number = numberInColumnPlus[iColumn];
//*** modify so sees if can go in
//see if it can go in at end
int * nextColumn = nextColumn_.array();
int * lastColumn = lastColumn_.array();
if (lengthAreaR_-startR[maximumColumnsExtra_]<number+1) {
//compression
int jColumn = nextColumn[maximumColumnsExtra_];
CoinBigIndex put = 0;
while ( jColumn != maximumColumnsExtra_ ) {
//move
CoinBigIndex get;
CoinBigIndex getEnd;
get = startR[jColumn];
getEnd = get + numberInColumnPlus[jColumn];
startR[jColumn] = put;
CoinBigIndex i;
for ( i = get; i < getEnd; i++ ) {
indexRowR[put] = indexRowR[i];
elementR[put] = elementR[i];
put++;
}
jColumn = nextColumn[jColumn];
}
numberCompressions_++;
startR[maximumColumnsExtra_]=put;
}
// Still may not be room (as iColumn was still in)
if (lengthAreaR_-startR[maximumColumnsExtra_]<number+1)
return false;
int next = nextColumn[iColumn];
int last = lastColumn[iColumn];
//out
nextColumn[last] = next;
lastColumn[next] = last;
CoinBigIndex put = startR[maximumColumnsExtra_];
//in at end
last = lastColumn[maximumColumnsExtra_];
nextColumn[last] = iColumn;
lastColumn[maximumColumnsExtra_] = iColumn;
lastColumn[iColumn] = last;
nextColumn[iColumn] = maximumColumnsExtra_;
//move
CoinBigIndex get = startR[iColumn];
startR[iColumn] = put;
int i = 0;
for (i=0 ; i < number; i ++ ) {
elementR[put]= elementR[get];
indexRowR[put++] = indexRowR[get++];
}
//insert
elementR[put]=value;
indexRowR[put++]=iRow;
numberInColumnPlus[iColumn]++;
//add 4 for luck
startR[maximumColumnsExtra_] = CoinMin((CoinBigIndex) (put + 4) ,lengthAreaR_);
return true;
}
/* getColumnSpaceIterate. Gets space for one extra U element in Column
may have to do compression (returns true)
also moves existing vector.
Returns -1 if no memory or where element was put
Used by replaceRow (turns off R version) */
CoinBigIndex
CoinFactorization::getColumnSpaceIterate ( int iColumn, double value,
int iRow)
{
if (numberInColumnPlus_.array()) {
numberInColumnPlus_.conditionalDelete();
}
int * numberInRow = numberInRow_.array();
int * numberInColumn = numberInColumn_.array();
int * nextColumn = nextColumn_.array();
int * lastColumn = lastColumn_.array();
int number = numberInColumn[iColumn];
int iNext=nextColumn[iColumn];
CoinBigIndex * startColumnU = startColumnU_.array();
CoinBigIndex * startRowU = startRowU_.array();
CoinBigIndex space = startColumnU[iNext]-startColumnU[iColumn];
CoinBigIndex put;
CoinBigIndex *convertRowToColumnU = convertRowToColumnU_.array();
int * indexColumnU = indexColumnU_.array();
double *elementU = elementU_.array();
int * indexRowU = indexRowU_.array();
if ( space < number + 1 ) {
//see if it can go in at end
if (lengthAreaU_-startColumnU[maximumColumnsExtra_]<number+1) {
//compression
int jColumn = nextColumn[maximumColumnsExtra_];
CoinBigIndex put = 0;
while ( jColumn != maximumColumnsExtra_ ) {
//move
CoinBigIndex get;
CoinBigIndex getEnd;
get = startColumnU[jColumn];
getEnd = get + numberInColumn[jColumn];
startColumnU[jColumn] = put;
CoinBigIndex i;
for ( i = get; i < getEnd; i++ ) {
double value = elementU[i];
if (value) {
indexRowU[put] = indexRowU[i];
elementU[put] = value;
put++;
} else {
numberInColumn[jColumn]--;
}
}
jColumn = nextColumn[jColumn];
}
numberCompressions_++;
startColumnU[maximumColumnsExtra_]=put;
//space for cross reference
CoinBigIndex *convertRowToColumn = convertRowToColumnU_.array();
CoinBigIndex j = 0;
CoinBigIndex *startRow = startRowU_.array();
int iRow;
for ( iRow = 0; iRow < numberRowsExtra_; iRow++ ) {
startRow[iRow] = j;
j += numberInRow[iRow];
}
factorElements_=j;
CoinZeroN ( numberInRow, numberRowsExtra_ );
int i;
for ( i = 0; i < numberRowsExtra_; i++ ) {
CoinBigIndex start = startColumnU[i];
CoinBigIndex end = start + numberInColumn[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;
}
}
}
// Still may not be room (as iColumn was still in)
if (lengthAreaU_-startColumnU[maximumColumnsExtra_]>=number+1) {
int next = nextColumn[iColumn];
int last = lastColumn[iColumn];
//out
nextColumn[last] = next;
lastColumn[next] = last;
put = startColumnU[maximumColumnsExtra_];
//in at end
last = lastColumn[maximumColumnsExtra_];
nextColumn[last] = iColumn;
lastColumn[maximumColumnsExtra_] = iColumn;
lastColumn[iColumn] = last;
nextColumn[iColumn] = maximumColumnsExtra_;
//move
CoinBigIndex get = startColumnU[iColumn];
startColumnU[iColumn] = put;
int i = 0;
for (i=0 ; i < number; i ++ ) {
double value = elementU[get];
int iRow=indexRowU[get++];
if (value) {
elementU[put]= value;
int n=numberInRow[iRow];
CoinBigIndex start = startRowU[iRow];
CoinBigIndex j;
for (j=start;j<start+n;j++) {
if (indexColumnU[j]==iColumn) {
convertRowToColumnU[j]=put;
break;
}
}
assert (j<start+n);
indexRowU[put++] = iRow;
} else {
assert (!numberInRow[iRow]);
numberInColumn[iColumn]--;
}
}
//insert
int n=numberInRow[iRow];
CoinBigIndex start = startRowU[iRow];
CoinBigIndex j;
for (j=start;j<start+n;j++) {
if (indexColumnU[j]==iColumn) {
convertRowToColumnU[j]=put;
break;
}
}
assert (j<start+n);
elementU[put]=value;
indexRowU[put]=iRow;
numberInColumn[iColumn]++;
//add 4 for luck
startColumnU[maximumColumnsExtra_] = CoinMin((CoinBigIndex) (put + 4) ,lengthAreaU_);
} else {
// no room
put=-1;
}
} else {
// just slot in
put=startColumnU[iColumn]+numberInColumn[iColumn];
int n=numberInRow[iRow];
CoinBigIndex start = startRowU[iRow];
CoinBigIndex j;
for (j=start;j<start+n;j++) {
if (indexColumnU[j]==iColumn) {
convertRowToColumnU[j]=put;
break;
}
}
assert (j<start+n);
elementU[put]=value;
indexRowU[put]=iRow;
numberInColumn[iColumn]++;
}
return put;
}
syntax highlighted by Code2HTML, v. 0.9.1