limp-cbc-0.3.2.0: cbits/coin/CoinOslFactorization.cpp
/* $Id: CoinOslFactorization.cpp 1585 2013-04-06 20:42:02Z stefan $ */
// Copyright (C) 1987, 2009, International Business Machines
// Corporation and others. All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).
#include "CoinUtilsConfig.h"
#include <cassert>
#include "CoinPragma.hpp"
#include "CoinOslFactorization.hpp"
#include "CoinOslC.h"
#include "CoinIndexedVector.hpp"
#include "CoinHelperFunctions.hpp"
#include "CoinPackedMatrix.hpp"
#include "CoinTypes.hpp"
#include "CoinFinite.hpp"
#include <stdio.h>
static void c_ekksmem(EKKfactinfo *fact,int numberRows,int maximumPivots);
static void c_ekksmem_copy(EKKfactinfo *fact,const EKKfactinfo * rhsFact);
static void c_ekksmem_delete(EKKfactinfo *fact);
//:class CoinOslFactorization. Deals with Factorization and Updates
// CoinOslFactorization. Constructor
CoinOslFactorization::CoinOslFactorization ( )
: CoinOtherFactorization()
{
gutsOfInitialize();
}
/// Copy constructor
CoinOslFactorization::CoinOslFactorization ( const CoinOslFactorization &other)
: CoinOtherFactorization(other)
{
gutsOfInitialize();
gutsOfCopy(other);
}
// Clone
CoinOtherFactorization *
CoinOslFactorization::clone() const
{
return new CoinOslFactorization(*this);
}
/// The real work of constructors etc
void CoinOslFactorization::gutsOfDestructor(bool clearFact)
{
delete [] elements_;
delete [] pivotRow_;
delete [] workArea_;
elements_ = NULL;
pivotRow_ = NULL;
workArea_ = NULL;
numberRows_ = 0;
numberColumns_ = 0;
numberGoodU_ = 0;
status_ = -1;
maximumRows_=0;
maximumSpace_=0;
solveMode_=0;
if (clearFact)
c_ekksmem_delete(&factInfo_);
}
void CoinOslFactorization::gutsOfInitialize(bool zapFact)
{
pivotTolerance_ = 1.0e-1;
zeroTolerance_ = 1.0e-13;
#ifndef COIN_FAST_CODE
slackValue_ = -1.0;
#endif
maximumPivots_=200;
relaxCheck_=1.0;
numberRows_ = 0;
numberColumns_ = 0;
numberGoodU_ = 0;
status_ = -1;
numberPivots_ = 0;
maximumRows_=0;
maximumSpace_=0;
elements_ = NULL;
pivotRow_ = NULL;
workArea_ = NULL;
solveMode_=0;
if (zapFact) {
memset(&factInfo_,0,sizeof(factInfo_));
factInfo_.maxinv=100;
factInfo_.drtpiv=1.0e-10;
factInfo_.zeroTolerance=1.0e-12;
factInfo_.zpivlu=0.1;
factInfo_.areaFactor=1.0;
factInfo_.nbfinv=100;
}
}
// ~CoinOslFactorization. Destructor
CoinOslFactorization::~CoinOslFactorization ( )
{
gutsOfDestructor();
}
// =
CoinOslFactorization & CoinOslFactorization::operator = ( const CoinOslFactorization & other ) {
if (this != &other) {
bool noGood = factInfo_.nrowmx!=other.factInfo_.nrowmx&&
factInfo_.eta_size!=other.factInfo_.eta_size;
gutsOfDestructor(noGood);
gutsOfInitialize(noGood);
gutsOfCopy(other);
}
return *this;
}
#define WORK_MULT 2
void CoinOslFactorization::gutsOfCopy(const CoinOslFactorization &other)
{
pivotTolerance_ = other.pivotTolerance_;
zeroTolerance_ = other.zeroTolerance_;
#ifndef COIN_FAST_CODE
slackValue_ = other.slackValue_;
#endif
relaxCheck_ = other.relaxCheck_;
numberRows_ = other.numberRows_;
numberColumns_ = other.numberColumns_;
maximumRows_ = other.maximumRows_;
maximumSpace_ = other.maximumSpace_;
solveMode_ = other.solveMode_;
numberGoodU_ = other.numberGoodU_;
maximumPivots_ = other.maximumPivots_;
numberPivots_ = other.numberPivots_;
factorElements_ = other.factorElements_;
status_ = other.status_;
elements_ = NULL;
pivotRow_ = NULL;
workArea_ = NULL;
c_ekksmem_copy(&factInfo_,&other.factInfo_);
}
// getAreas. Gets space for a factorization
//called by constructors
void
CoinOslFactorization::getAreas ( int numberOfRows,
int numberOfColumns,
CoinBigIndex maximumL,
CoinBigIndex maximumU )
{
numberRows_ = numberOfRows;
numberColumns_ = numberOfColumns;
CoinBigIndex size = static_cast<CoinBigIndex>(factInfo_.areaFactor*
(maximumL+maximumU));
factInfo_.zeroTolerance=zeroTolerance_;
// If wildly out redo
if (maximumRows_>numberRows_+1000) {
maximumRows_=0;
maximumSpace_=0;
factInfo_.last_eta_size=0;
}
if (size>maximumSpace_) {
//delete [] elements_;
//elements_ = new CoinFactorizationDouble [size];
maximumSpace_ = size;
}
factInfo_.lastEtaCount = factInfo_.nnentu+factInfo_.nnentl;
int oldnnetas=factInfo_.last_eta_size;
// If we are going to increase then be on safe side
if (size>oldnnetas)
size = static_cast<int>(1.1*size);
factInfo_.eta_size=CoinMax(size,oldnnetas);
//printf("clp size %d, old %d now %d - iteration %d - last count %d - rows %d,%d,%d\n",
// size,oldnnetas,factInfo_.eta_size,factInfo_.iterno,factInfo_.lastEtaCount,
//numberRows_,factInfo_.nrowmx,factInfo_.nrow);
//if (!factInfo_.iterno) {
//printf("here\n");
//}
/** Get solve mode e.g. 0 C++ code, 1 Lapack, 2 choose
If 4 set then values pass
if 8 set then has iterated
*/
solveMode_ &= 4+8; // clear bottom bits
factInfo_.ifvsol= ((solveMode_&4)!=0) ? 1 : 0;
if ((solveMode_&8)!=0) {
factInfo_.ifvsol=0;
factInfo_.invok=1;
} else {
factInfo_.iter0=factInfo_.iterno;
factInfo_.invok=-1;
factInfo_.if_sparse_update=0;
}
#if 0
if (!factInfo_.if_sparse_update &&
factInfo_.iterno>factInfo_.iter0 &&
numberRows_>=C_EKK_GO_SPARSE) {
printf("count %d rows %d etasize %d\n",
factInfo_.lastEtaCount,factInfo_.nrow,factInfo_.eta_size);
}
#endif
if (!factInfo_.if_sparse_update &&
factInfo_.iterno>factInfo_.iter0 &&
numberRows_>=C_EKK_GO_SPARSE &&
(factInfo_.lastEtaCount>>2)<factInfo_.nrow&&
!factInfo_.switch_off_sparse_update) {
#if PRINT_DEBUG
printf("**** Switching on sparse update - etacount\n");
#endif
/* I suspect this can go into c_ekksslvf;
* if c_ekkshff decides to switch sparse_update off,
* then it problably always switches it off (?)
*/
factInfo_.if_sparse_update=2;
}
c_ekksmem(&factInfo_,numberRows_,maximumPivots_);
if (numberRows_>maximumRows_) {
maximumRows_ = numberRows_;
//delete [] pivotRow_;
//delete [] workArea_;
//pivotRow_ = new int [2*maximumRows_+maximumPivots_];
//workArea_ = new CoinFactorizationDouble [maximumRows_*WORK_MULT];
}
}
// preProcess.
void
CoinOslFactorization::preProcess ()
{
factInfo_.zpivlu=pivotTolerance_;
// Go to Fortran
int * hcoli=factInfo_.xecadr+1;
int * indexRowU = factInfo_.xeradr+1;
CoinBigIndex * startColumnU=factInfo_.xcsadr+1;
for (int i=0;i<numberRows_;i++) {
int start = startColumnU[i];
startColumnU[i]++; // to Fortran
for (int j=start;j<startColumnU[i+1];j++) {
indexRowU[j]++; // to Fortran
hcoli[j]=i+1; // to Fortran
}
}
startColumnU[numberRows_]++; // to Fortran
/* can do in column order - no zeros or duplicates */
#ifndef NDEBUG
int ninbas =
#endif
c_ekkslcf(&factInfo_);
assert (ninbas>0);
}
//Does factorization
int
CoinOslFactorization::factor ( )
{
/* Uwe's factorization (sort of) */
int irtcod = c_ekklfct(&factInfo_);
/* Check return code */
/* 0 - Fine , 1 - Backtrack, 2 - Singularities on initial, 3-Fatal */
/* now 5-Need more memory */
status_= 0;
if (factInfo_.eta_size>factInfo_.last_eta_size) {
factInfo_.areaFactor *= factInfo_.eta_size;
factInfo_.areaFactor /= factInfo_.last_eta_size;
#ifdef CLP_INVESTIGATE
printf("areaFactor increased to %g\n",factInfo_.areaFactor);
#endif
}
if (irtcod==5) {
status_=-99;
assert (factInfo_.eta_size>factInfo_.last_eta_size) ;
#ifdef CLP_INVESTIGATE
printf("need more memory\n");
#endif
} else if (irtcod) {
status_=-1;
//printf("singular %d\n",irtcod);
}
return status_;
}
// Makes a non-singular basis by replacing variables
void
CoinOslFactorization::makeNonSingular(int * sequence, int numberColumns)
{
const EKKHlink *rlink = factInfo_.kp1adr;
const EKKHlink *clink = factInfo_.kp2adr;
int nextRow=0;
//int * mark = reinterpret_cast<int *>(factInfo_.kw1adr);
//int nr=0;
//int nc=0;
#if 0
for (int i=0;i<numberRows_;i++) {
//mark[i]=-1;
if (rlink[i].pre>=0||rlink[i].pre==-(numberRows_+1)) {
nr++;
printf("%d rl %d cl %d\n",i,rlink[i].pre,clink[i].pre);
}
if (clink[i].pre>=0||clink[i].pre==-(numberRows_+1)) {
nc++;
printf("%d rl %d cl %d\n",i,rlink[i].pre,clink[i].pre);
}
}
#endif
//printf("nr %d nc %d\n",nr,nc);
#ifndef NDEBUG
bool goodPass=true;
#endif
int numberDone=0;
for (int i=0;i<numberRows_;i++) {
int cRow =(-clink[i].pre)-1;
if (cRow==numberRows_||cRow<0) {
// throw out
for (;nextRow<numberRows_;nextRow++) {
int rRow =(-rlink[nextRow].pre)-1;
if (rRow==numberRows_||rRow<0)
break;
}
if (nextRow<numberRows_) {
sequence[i]=nextRow+numberColumns;
nextRow++;
numberDone++;
} else {
#ifndef NDEBUG
goodPass=false;
#endif
assert(numberDone);
//printf("BAD singular at row %d\n",i);
break;
}
}
}
#ifndef NDEBUG
if (goodPass) {
for (;nextRow<numberRows_;nextRow++) {
int rRow =(-rlink[nextRow].pre)-1;
assert (!(rRow==numberRows_||rRow<0));
}
}
#endif
}
// Does post processing on valid factorization - putting variables on correct rows
void
CoinOslFactorization::postProcess(const int * sequence, int * pivotVariable)
{
factInfo_.iterin=factInfo_.iterno;
factInfo_.npivots=0;
numberPivots_=0;
const int * permute3 = factInfo_.mpermu+1;
assert (permute3==reinterpret_cast<const int *>
(factInfo_.kadrpm+numberRows_+1));
// this is ridiculous - must be better way
int * permute2 = reinterpret_cast<int *>(factInfo_.kw1adr);
const int * permute = reinterpret_cast<const int *>(factInfo_.kp2adr);
for (int i=0;i<numberRows_;i++) {
permute2[permute[i]-1]=i;
}
for (int i=0;i<numberRows_;i++) {
// the row is i
// the column is whatever matches k3[i] in k1
int look=permute3[i]-1;
int j=permute2[look];
int k = sequence[j];
pivotVariable[i]=k;
}
#ifdef CLP_REUSE_ETAS
int * start = factInfo_.xcsadr+1;
int * putSeq = factInfo_.xrsadr+2*factInfo_.nrowmx+2;
int * position = putSeq+factInfo_.maxinv;
int * putStart = position+factInfo_.maxinv;
memcpy(putStart,start,numberRows_*sizeof(int));
int iLast=start[numberRows_-1];
putStart[numberRows_]=iLast+factInfo_.xeradr[iLast]+1;
factInfo_.save_nnentu=factInfo_.nnentu;
#endif
#ifndef NDEBUG
{
int lstart=numberRows_+factInfo_.maxinv+5;
int ndo = factInfo_.xnetal-lstart;
double * dluval=factInfo_.xeeadr;
int * mcstrt = factInfo_.xcsadr+lstart;
if (ndo)
assert (dluval[mcstrt[ndo]+1]<1.0e50);
}
#endif
}
/* Replaces one Column to basis,
returns 0=OK, 1=Probably OK, 2=singular, 3=no room
If checkBeforeModifying is true will do all accuracy checks
before modifying factorization. Whether to set this depends on
speed considerations. You could just do this on first iteration
after factorization and thereafter re-factorize
partial update already in U */
int
CoinOslFactorization::replaceColumn ( CoinIndexedVector * regionSparse,
int pivotRow,
double pivotCheck ,
bool /*checkBeforeModifying*/,
double acceptablePivot)
{
if (numberPivots_+1==maximumPivots_)
return 3;
int *regionIndex = regionSparse->getIndices ( );
double *region = regionSparse->denseVector ( );
int orig_nincol=0;
double saveTolerance = factInfo_.drtpiv;
factInfo_.drtpiv=acceptablePivot;
int returnCode=c_ekketsj(&factInfo_,region-1,
regionIndex,
pivotCheck,orig_nincol,
numberPivots_,&factInfo_.nuspike,
pivotRow+1,
reinterpret_cast<int *>(factInfo_.kw1adr));
factInfo_.drtpiv=saveTolerance;
if (returnCode!=2)
numberPivots_++;
#ifndef NDEBUG
{
int lstart=numberRows_+factInfo_.maxinv+5;
int ndo = factInfo_.xnetal-lstart;
double * dluval=factInfo_.xeeadr;
int * mcstrt = factInfo_.xcsadr+lstart;
if (ndo)
assert (dluval[mcstrt[ndo]+1]<1.0e50);
}
#endif
return returnCode;
}
/* This version has same effect as above with FTUpdate==false
so number returned is always >=0 */
int
CoinOslFactorization::updateColumn ( CoinIndexedVector * regionSparse,
CoinIndexedVector * regionSparse2,
bool /*noPermute*/) const
{
#ifndef NDEBUG
{
int lstart=numberRows_+factInfo_.maxinv+5;
int ndo = factInfo_.xnetal-lstart;
double * dluval=factInfo_.xeeadr;
int * mcstrt = factInfo_.xcsadr+lstart;
if (ndo)
assert (dluval[mcstrt[ndo]+1]<1.0e50);
}
#endif
assert (numberRows_==numberColumns_);
double *region2 = regionSparse2->denseVector ( );
int *regionIndex2 = regionSparse2->getIndices ( );
int numberNonZero = regionSparse2->getNumElements ( );
double *region = regionSparse->denseVector ( );
//const int * permuteIn = factInfo_.mpermu+1;
// Stuff is put one up so won't get illegal read
assert (!region[numberRows_]);
assert (!regionSparse2->packedMode());
#if 0
int first=numberRows_;
for (int j=0;j<numberNonZero;j++) {
int jRow = regionIndex2[j];
int iRow = permuteIn[jRow];
region[iRow]=region2[jRow];
first=CoinMin(first,iRow);
region2[jRow]=0.0;
}
#endif
numberNonZero=c_ekkftrn(&factInfo_,
region2-1,region,regionIndex2,numberNonZero);
regionSparse2->setNumElements(numberNonZero);
return 0;
}
/* Updates one column (FTRAN) from regionSparse2
Tries to do FT update
number returned is negative if no room
regionSparse starts as zero and is zero at end.
Note - if regionSparse2 packed on input - will be packed on output
*/
int
CoinOslFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
CoinIndexedVector * regionSparse2,
bool /*noPermute*/)
{
assert (numberRows_==numberColumns_);
double *region2 = regionSparse2->denseVector ( );
int *regionIndex2 = regionSparse2->getIndices ( );
int numberNonZero = regionSparse2->getNumElements ( );
assert (regionSparse2->packedMode());
// packed mode
//int numberNonInOriginal=numberNonZero;
//double *dpermu = factInfo_.kadrpm;
// Use region instead of dpermu
double * save =factInfo_.kadrpm;
factInfo_.kadrpm=regionSparse->denseVector()-1;
int nuspike=c_ekkftrn_ft(&factInfo_, region2,regionIndex2,
&numberNonZero);
factInfo_.kadrpm=save;
regionSparse2->setNumElements(numberNonZero);
//regionSparse2->print();
factInfo_.nuspike=nuspike;
return nuspike;
}
int
CoinOslFactorization::updateTwoColumnsFT(CoinIndexedVector * regionSparse1,
CoinIndexedVector * regionSparse2,
CoinIndexedVector * regionSparse3,
bool /*noPermute*/)
{
#if 1
// probably best to merge on a LU part by part
// but can try full merge
double *region2 = regionSparse2->denseVector ( );
int *regionIndex2 = regionSparse2->getIndices ( );
int numberNonZero2 = regionSparse2->getNumElements ( );
assert (regionSparse2->packedMode());
assert (numberRows_==numberColumns_);
double *region3 = regionSparse3->denseVector ( );
int *regionIndex3 = regionSparse3->getIndices ( );
int numberNonZero3 = regionSparse3->getNumElements ( );
double *region = regionSparse1->denseVector ( );
// Stuff is put one up so won't get illegal read
assert (!region[numberRows_]);
assert (!regionSparse3->packedMode());
// packed mode
//double *dpermu = factInfo_.kadrpm;
#if 0
factInfo_.nuspike=c_ekkftrn_ft(&factInfo_, region2,regionIndex2,
&numberNonZero2);
numberNonZero3=c_ekkftrn(&factInfo_,
region3-1,region,regionIndex3,numberNonZero3);
#else
c_ekkftrn2(&factInfo_,region3-1,region,regionIndex3,&numberNonZero3,
region2,regionIndex2,&numberNonZero2);
#endif
regionSparse2->setNumElements(numberNonZero2);
regionSparse3->setNumElements(numberNonZero3);
return factInfo_.nuspike;
#else
// probably best to merge on a LU part by part
// but can try full merge
int returnCode= updateColumnFT(regionSparse1,
regionSparse2);
updateColumn(regionSparse1,
regionSparse3,
noPermute);
return returnCode;
#endif
}
/* Updates one column (BTRAN) from regionSparse2
regionSparse starts as zero and is zero at end
Note - if regionSparse2 packed on input - will be packed on output
*/
int
CoinOslFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
CoinIndexedVector * regionSparse2) const
{
assert (numberRows_==numberColumns_);
double *region2 = regionSparse2->denseVector ( );
int *regionIndex2 = regionSparse2->getIndices ( );
int numberNonZero = regionSparse2->getNumElements ( );
//double *region = regionSparse->denseVector ( );
/*int *regionIndex = regionSparse->getIndices ( );*/
const int * permuteIn = factInfo_.mpermu+1;
factInfo_.packedMode = regionSparse2->packedMode() ? 1 : 0;
// Use region instead of dpermu
double * save =factInfo_.kadrpm;
factInfo_.kadrpm=regionSparse->denseVector()-1;
// use internal one for now (address is one off)
double * region = factInfo_.kadrpm;
if (numberNonZero<2) {
if (numberNonZero) {
int ipivrw=regionIndex2[0];
if (factInfo_.packedMode) {
double value=region2[0];
region2[0]=0.0;
region2[ipivrw]=value;
}
numberNonZero=c_ekkbtrn_ipivrw(&factInfo_, region2-1,
regionIndex2-1,ipivrw+1,
reinterpret_cast<int *>(factInfo_.kp1adr));
}
} else {
#ifndef NDEBUG
{
int *mcstrt = factInfo_.xcsadr;
int * hpivco_new=factInfo_.kcpadr+1;
int nrow=factInfo_.nrow;
int i;
int ipiv = hpivco_new[0];
int last = mcstrt[ipiv];
for (i=0;i<nrow-1;i++) {
ipiv=hpivco_new[ipiv];
assert (mcstrt[ipiv]>last);
last=mcstrt[ipiv];
}
}
#endif
int iSmallest = COIN_INT_MAX;
int iPiv=0;
const int *mcstrt = factInfo_.xcsadr;
// permute and save where nonzeros are
if (!factInfo_.packedMode) {
if ((numberRows_<200||(numberNonZero<<4)>numberRows_)) {
for (int j=0;j<numberNonZero;j++) {
int jRow = regionIndex2[j];
int iRow = permuteIn[jRow];
regionIndex2[j]=iRow;
region[iRow]=region2[jRow];
region2[jRow]=0.0;
}
} else {
for (int j=0;j<numberNonZero;j++) {
int jRow = regionIndex2[j];
int iRow = permuteIn[jRow];
regionIndex2[j]=iRow;
region[iRow]=region2[jRow];
if (mcstrt[iRow]<iSmallest) {
iPiv=iRow;
iSmallest=mcstrt[iRow];
}
region2[jRow]=0.0;
}
}
} else {
for (int j=0;j<numberNonZero;j++) {
int jRow = regionIndex2[j];
int iRow = permuteIn[jRow];
regionIndex2[j]=iRow;
region[iRow]=region2[j];
region2[j]=0.0;
}
}
assert (iPiv>=0);
numberNonZero=c_ekkbtrn(&factInfo_, region2-1,regionIndex2-1,iPiv);
}
factInfo_.kadrpm=save;
factInfo_.packedMode=0;
regionSparse2->setNumElements(numberNonZero);
return 0;
}
// Number of entries in each row
int *
CoinOslFactorization::numberInRow() const
{ return reinterpret_cast<int *> (factInfo_.xrnadr+1);}
// Number of entries in each column
int *
CoinOslFactorization::numberInColumn() const
{ return reinterpret_cast<int *> (factInfo_.xcnadr+1);}
// Returns array to put basis starts in
CoinBigIndex *
CoinOslFactorization::starts() const
{ return reinterpret_cast<CoinBigIndex *> (factInfo_.xcsadr+1);}
// Returns array to put basis elements in
CoinFactorizationDouble *
CoinOslFactorization::elements() const
{ return factInfo_.xeeadr+1;}
// Returns pivot row
int *
CoinOslFactorization::pivotRow() const
{ return factInfo_.krpadr+1;}
// Returns work area
CoinFactorizationDouble *
CoinOslFactorization::workArea() const
{ return factInfo_.kw1adr;}
// Returns int work area
int *
CoinOslFactorization::intWorkArea() const
{ return reinterpret_cast<int *> (factInfo_.kw1adr);}
// Returns permute back
int *
CoinOslFactorization::permuteBack() const
{ return factInfo_.kcpadr+1;}
// Returns array to put basis indices in
int *
CoinOslFactorization::indices() const
{ return factInfo_.xeradr+1;}
// Returns true if wants tableauColumn in replaceColumn
bool
CoinOslFactorization::wantsTableauColumn() const
{ return false;}
/* Useful information for factorization
0 - iteration number
whereFrom is 0 for factorize and 1 for replaceColumn
*/
#ifdef CLP_REUSE_ETAS
void
CoinOslFactorization::setUsefulInformation(const int * info,int whereFrom)
{
factInfo_.iterno=info[0];
if (whereFrom) {
factInfo_.reintro=-1;
if( factInfo_.first_dense>=factInfo_.last_dense) {
int * putSeq = factInfo_.xrsadr+2*factInfo_.nrowmx+2;
int * position = putSeq+factInfo_.maxinv;
//int * putStart = position+factInfo_.maxinv;
int iSequence=info[1];
if (whereFrom==1) {
putSeq[factInfo_.npivots]=iSequence;
} else {
int i;
for (i=factInfo_.npivots-1;i>=0;i--) {
if (putSeq[i]==iSequence)
break;
}
if (i>=0) {
factInfo_.reintro=position[i];
} else {
factInfo_.reintro=-1;
}
factInfo_.nnentu=factInfo_.save_nnentu;
}
}
}
}
#else
void
CoinOslFactorization::setUsefulInformation(const int * info,int /*whereFrom*/)
{ factInfo_.iterno=info[0]; }
#endif
// Get rid of all memory
void
CoinOslFactorization::clearArrays()
{
factInfo_.nR_etas=0;
factInfo_.nnentu=0;
factInfo_.nnentl=0;
maximumRows_=0;
maximumSpace_=0;
factInfo_.last_eta_size=0;
gutsOfDestructor(false);
}
void
CoinOslFactorization::maximumPivots ( int value )
{
maximumPivots_ = value;
}
#define CLP_FILL 15
/*#undef NDEBUG*/
//#define CLP_DEBUG_MALLOC 1000000
#if CLP_DEBUG_MALLOC
static int malloc_number=0;
static int malloc_check=-1;
static int malloc_counts_on=0;
struct malloc_struct {
void * previous;
void * next;
int size;
int when;
int type;
};
static double malloc_times=0.0;
static double malloc_total=0.0;
static double malloc_current=0.0;
static double malloc_max=0.0;
static int malloc_amount[]={0,32,128,256,1024,4096,16384,65536,262144,
2000000000};
static int malloc_n=10;
double malloc_counts[10]={0,0,0,0,0,0,0,0,0,0};
typedef struct malloc_struct malloc_struct;
static malloc_struct startM = {0,0,0,0};
static malloc_struct endM = {0,0,0,0};
static int extra=4;
void clp_memory(int type)
{
if (type==0) {
/* switch on */
malloc_counts_on=1;
startM.next=&endM;
endM.previous=&startM;
} else {
/* summary */
double average = malloc_total/malloc_times;
int i;
malloc_struct * previous = (malloc_struct *) endM.previous;
printf("count %g bytes %g - average %g\n",malloc_times,malloc_total,average);
printf("current bytes %g - maximum %g\n",malloc_current,malloc_max);
for ( i=0;i<malloc_n;i++)
printf("%g ",malloc_counts[i]);
printf("\n");
malloc_counts_on=0;
if (previous->previous!=&startM) {
int n=0;
printf("Allocated blocks\n");
while (previous->previous!=&startM) {
printf("(%d at %d) ",previous->size,previous->when);
n++;
if ((n%5)==0)
printf("\n");
previous = (malloc_struct *) previous->previous;
}
printf("\n - total %d\n",n);
}
}
malloc_number=0;
malloc_times=0.0;
malloc_total=0.0;
malloc_current=0.0;
malloc_max=0.0;
memset(malloc_counts,0,sizeof(malloc_counts));
}
static void clp_adjust(void * temp,int size,int type)
{
malloc_struct * itemp = (malloc_struct *) temp;
malloc_struct * first = (malloc_struct *) startM.next;
int i;
startM.next=temp;
first->previous=temp;
itemp->previous=&startM;
itemp->next=first;
malloc_number++;
if (malloc_number==malloc_check) {
printf("Allocation of %d bytes at %d (type %d) not freed\n",
size,malloc_number,type);
}
itemp->when=malloc_number;
itemp->size=size;
itemp->type=type;
malloc_times ++;
malloc_total += size;
malloc_current += size;
malloc_max=CoinMax(malloc_max,malloc_current);
for (i=0;i<malloc_n;i++) {
if ((int) size<=malloc_amount[i]) {
malloc_counts[i]++;
break;
}
}
}
#endif
/* covers */
double * clp_double(int number_entries)
{
#if CLP_DEBUG_MALLOC==0
return reinterpret_cast<double *>( malloc(number_entries*sizeof(double)));
#else
double * temp = reinterpret_cast<double *>( malloc((number_entries+extra)*sizeof(double)));
clp_adjust(temp,number_entries*sizeof(double),1);
#if CLP_DEBUG_MALLOC>1
if (number_entries*sizeof(double)>=CLP_DEBUG_MALLOC)
printf("WWW %x malloced by double %d - size %d\n",
temp+extra,malloc_number,number_entries);
#endif
return temp+extra;
#endif
}
int * clp_int(int number_entries)
{
#if CLP_DEBUG_MALLOC==0
return reinterpret_cast<int *>( malloc(number_entries*sizeof(int)));
#else
double * temp = reinterpret_cast<double *>( malloc(((number_entries+1)/2+extra)*sizeof(double)));
clp_adjust(temp,number_entries*sizeof(int),2);
#if CLP_DEBUG_MALLOC>1
if (number_entries*sizeof(int)>=CLP_DEBUG_MALLOC)
printf("WWW %x malloced by int %d - size %d\n",
temp+extra,malloc_number,number_entries);
#endif
return reinterpret_cast<int *>( (temp+extra));
#endif
}
void * clp_malloc(int number_entries)
{
#if CLP_DEBUG_MALLOC==0
return malloc(number_entries);
#else
double * temp = reinterpret_cast<double *>( malloc(number_entries+extra*sizeof(double)));
clp_adjust(temp,number_entries,0);
#if CLP_DEBUG_MALLOC>1
if (number_entries>=CLP_DEBUG_MALLOC)
printf("WWW %x malloced by void %d - size %d\n",
temp+extra,malloc_number,number_entries);
#endif
return (void *) (temp+extra);
#endif
}
void clp_free(void * oldArray)
{
#if CLP_DEBUG_MALLOC==0
free(oldArray);
#else
if (oldArray) {
double * temp = (reinterpret_cast<double *>( oldArray)-extra);
malloc_struct * itemp = (malloc_struct *) temp;
malloc_struct * next = (malloc_struct *) itemp->next;
malloc_struct * previous = (malloc_struct *) itemp->previous;
previous->next=next;
next->previous=previous;
malloc_current -= itemp->size;
#if CLP_DEBUG_MALLOC>1
if (itemp->size>=CLP_DEBUG_MALLOC)
printf("WWW %x freed by free %d - old length %d - type %d\n",
oldArray,itemp->when,itemp->size,itemp->type);
#endif
free(temp);
}
#endif
}
/*#define FIX_ADD 4*nrowmx+5
#define FIX_ADD2 4*nrowmx+5*/
#define FIX_ADD 1*nrowmx+5
#define FIX_ADD2 1*nrowmx+5
#define ALIGNMENT 32
static void * clp_align (void * memory)
{
if (sizeof(int)==sizeof(void *)&&ALIGNMENT) {
CoinInt64 k = reinterpret_cast<CoinInt64> (memory);
if ((k&(ALIGNMENT-1))!=0) {
k &= ~(ALIGNMENT-1);
k += ALIGNMENT;
memory = reinterpret_cast<void *> (k);
}
return memory;
} else {
return memory;
}
}
void clp_setup_pointers(EKKfactinfo * fact)
{
/* do extra stuff */
int nrow=fact->nrow;
int nrowmx=fact->nrowmx;
int maxinv=fact->maxinv;
fact->lstart = nrow + maxinv + 5;
/* this is the number of L transforms */
fact->xnetalval = fact->xnetal - fact->lstart;
fact->mpermu = (reinterpret_cast<int*> (fact->kadrpm+nrow))+1;
fact->bitArray = fact->krpadr + ( nrowmx+2);
fact->back = fact->kcpadr+2*nrow + maxinv + 4;
fact->hpivcoR = fact->kcpadr+nrow+3;
fact->nonzero = (reinterpret_cast<char *>( &fact->mpermu[nrow+1]))-1;
}
#ifndef NDEBUG
int ets_count=0;
int ets_check=-1;
//static int adjust_count=0;
//static int adjust_check=-1;
#endif
static void clp_adjust_pointers(EKKfactinfo * fact, int adjust)
{
#if 0 //ndef NDEBUG
adjust_count++;
if (adjust_check>=0&&adjust_count>=adjust_check) {
printf("trouble\n");
}
#endif
if (fact->trueStart) {
fact->kadrpm += adjust;
fact->krpadr += adjust;
fact->kcpadr += adjust;
fact->xrsadr += adjust;
fact->xcsadr += adjust;
fact->xrnadr += adjust;
fact->xcnadr += adjust;
}
if (fact->xeradr) {
fact->xeradr += adjust;
fact->xecadr += adjust;
fact->xeeadr += adjust;
}
}
/* deals with memory for complicated array
0 just do addresses
1 just get memory */
static double *
clp_alloc_memory(EKKfactinfo * fact,int type, int * length)
{
int nDouble=0;
int nInt=0;
int nrowmxp;
int ntot1;
int ntot2;
int ntot3;
int nrowmx;
int * tempI;
double * tempD;
nrowmx=fact->nrowmx;
nrowmxp = nrowmx + 2;
ntot1 = nrowmxp;
ntot2 = 3*nrowmx+5; /* space for three lists */
ntot3 = 2*nrowmx;
if ((ntot1<<1)<ntot2) {
ntot1=ntot2>>1;
}
ntot3=CoinMax(ntot3,ntot1);
/* Row work regions */
/* must be contiguous so allocate as one chunk */
/* may only need 2.5 */
/* now doing all at once - far too much - reduce later */
tempD=fact->kw1adr;
tempD+=nrowmxp;
tempD = reinterpret_cast<double *>( clp_align(tempD));
fact->kw2adr=tempD;
tempD+=nrowmxp;
tempD = reinterpret_cast<double *>( clp_align(tempD));
fact->kw3adr=tempD-1;
tempD+=nrowmxp;
tempD = reinterpret_cast<double *>( clp_align(tempD));
fact->kp1adr=reinterpret_cast<EKKHlink *>(tempD);
tempD+=nrowmxp;
tempD = reinterpret_cast<double *>( clp_align(tempD));
fact->kp2adr=reinterpret_cast<EKKHlink *>(tempD);
//tempD+=ntot3;
tempD+=nrowmxp;
tempD = reinterpret_cast<double *>( clp_align(tempD));
/*printf("zz %x %x\n",tempD,fact->kadrpm);*/
fact->kadrpm = tempD;
/* seems a lot */
tempD += ((6*nrowmx +8)*(sizeof(int))/sizeof(double));
/* integer arrays */
tempI = reinterpret_cast<int *>( tempD);
tempI = reinterpret_cast<int *>( clp_align(tempI));
fact->xrsadr = tempI;
#ifdef CLP_REUSE_ETAS
tempI +=( 3*(nrowmx+fact->maxinv+1));
#else
tempI +=( (nrowmx<<1)+fact->maxinv+1);
#endif
tempI = reinterpret_cast<int *>( clp_align(tempI));
fact->xcsadr = tempI;
#if 1 //def CLP_REUSE_ETAS
tempI += ( 2*nrowmx+8+2*fact->maxinv);
#else
tempI += ( 2*nrowmx+8+fact->maxinv);
#endif
tempI += FIX_ADD+FIX_ADD2;
tempI = reinterpret_cast<int *>( clp_align(tempI));
fact->xrnadr = tempI;
tempI += nrowmx;
tempI = reinterpret_cast<int *>( clp_align(tempI));
fact->xcnadr = tempI;
tempI += nrowmx;
tempI = reinterpret_cast<int *>( clp_align(tempI));
fact->krpadr = tempI;
tempI += ( nrowmx+1) +((nrowmx+33)>>5);
/*printf("zzz %x %x\n",tempI,fact->kcpadr);*/
tempI = reinterpret_cast<int *>( clp_align(tempI));
fact->kcpadr = tempI;
tempI += 3*nrowmx+8+fact->maxinv;
fact->R_etas_start = fact->xcsadr+nrowmx+fact->maxinv+4;
fact->R_etas_start += FIX_ADD;
nInt = static_cast<int>(tempI-(reinterpret_cast<int *>( fact->trueStart)));
nDouble = static_cast<int>(sizeof(int)*(nInt+1)/sizeof(double));
*length = nDouble;
/*printf("nDouble %d - type %d\n",nDouble,type);*/
nDouble += static_cast<int>((2*ALIGNMENT)/sizeof(double));
if (type) {
/*printf("%d allocated\n",nDouble);*/
tempD = reinterpret_cast<double *>( clp_double(nDouble));
#ifndef NDEBUG
memset(tempD,CLP_FILL,nDouble*sizeof(double));
#endif
}
return tempD;
}
static void c_ekksmem(EKKfactinfo *fact,int nrow,int maximumPivots)
{
/* space for invert */
int nnetas=fact->eta_size;
fact->nrow=nrow;
if (!(nnetas>fact->last_eta_size||(!fact->xe2adr&&fact->if_sparse_update)||
nrow>fact->nrowmx||maximumPivots>fact->maxinv))
return;
clp_adjust_pointers(fact, +1);
if (nrow>fact->nrowmx||maximumPivots>fact->maxinv) {
int length;
fact->nrowmx=CoinMax(nrow,fact->nrowmx);
fact->maxinv=CoinMax(maximumPivots,fact->maxinv);
clp_free(fact->trueStart);
fact->trueStart=0;
fact->kw1adr=0;
fact->trueStart=clp_alloc_memory(fact,1,&length);
fact->kw1adr=reinterpret_cast<double *>( clp_align(fact->trueStart));
clp_alloc_memory(fact,0,&length);
}
/*if (!fact->iterno) fact->eta_size+=1000000;*//* TEMP*/
if (nnetas>fact->last_eta_size||(!fact->xe2adr&&fact->if_sparse_update)) {
fact->last_eta_size = nnetas;
clp_free(reinterpret_cast<char *>(fact->xe2adr));
/* if malloc fails - we have lost memory - start again */
if (!fact->ndenuc &&fact->if_sparse_update) {
/* allow second copy of elements */
fact->xe2adr = clp_double(nnetas);
#ifndef NDEBUG
memset(fact->xe2adr,CLP_FILL,nnetas*sizeof(double));
#endif
if (!fact->xe2adr) {
fact->maxNNetas=fact->last_eta_size; /* dont allow any increase */
nnetas=fact->last_eta_size;
fact->eta_size=nnetas;
#ifdef PRINT_DEBUG
if (fact->if_sparse_update) {
printf("*** Sparse update off due to memory\n");
}
#endif
fact->if_sparse_update=0;
fact->switch_off_sparse_update=1;
}
} else {
fact->xe2adr = 0;
fact->if_sparse_update=0;
}
clp_free(fact->xeradr);
fact->xeradr= clp_int( nnetas);
#ifndef NDEBUG
memset(fact->xeradr,CLP_FILL,nnetas*sizeof(int));
#endif
if (!fact->xeradr) {
nnetas=0;
}
if (nnetas) {
clp_free(fact->xecadr);
fact->xecadr= clp_int( nnetas);
#ifndef NDEBUG
memset(fact->xecadr,CLP_FILL,nnetas*sizeof(int));
#endif
if (!fact->xecadr) {
nnetas=0;
}
}
if (nnetas) {
clp_free(fact->xeeadr);
fact->xeeadr= clp_double(nnetas);
#ifndef NDEBUG
memset(fact->xeeadr,CLP_FILL,nnetas*sizeof(double));
#endif
if (!fact->xeeadr) {
nnetas=0;
}
}
}
if (!nnetas) {
char msg[100];
sprintf(msg,"Unable to allocate factorization memory for %d elements",
nnetas);
throw(msg);
}
/*c_ekklplp->nnetas=nnetas;*/
fact->nnetas=nnetas;
clp_adjust_pointers(fact, -1);
}
static void c_ekksmem_copy(EKKfactinfo *fact,const EKKfactinfo * rhsFact)
{
/* space for invert */
int nrowmx=rhsFact->nrowmx,nnetas=rhsFact->nnetas;
int canReuseEtas= (fact->eta_size==rhsFact->eta_size) ? 1 : 0;
int canReuseArrays = (fact->nrowmx==rhsFact->nrowmx) ? 1 : 0;
clp_adjust_pointers(fact, +1);
clp_adjust_pointers(const_cast<EKKfactinfo *>(rhsFact), +1);
/*memset(fact,0,sizeof(EKKfactinfo));*/
/* copy scalars */
memcpy(&fact->drtpiv,&rhsFact->drtpiv,5*sizeof(double));
memcpy(&fact->nrow,&rhsFact->nrow,((&fact->maxNNetas-&fact->nrow)+1)*
sizeof(int));
if (nrowmx) {
int length;
int kCopyEnd,nCopyEnd,nCopyStart;
if (!canReuseEtas) {
clp_free(fact->xeradr);
clp_free(fact->xecadr);
clp_free(fact->xeeadr);
clp_free(fact->xe2adr);
fact->xeradr = 0;
fact->xecadr = 0;
fact->xeeadr = 0;
fact->xe2adr = 0;
}
if (!canReuseArrays) {
clp_free(fact->trueStart);
fact->trueStart=0;
fact->kw1adr=0;
fact->trueStart=clp_alloc_memory(fact,1,&length);
fact->kw1adr=reinterpret_cast<double *>( clp_align(fact->trueStart));
}
clp_alloc_memory(fact,0,&length);
nnetas=fact->eta_size;
assert (nnetas);
{
int n2 = rhsFact->nR_etas;
int n3 = n2 ? rhsFact->R_etas_start[1+n2]: 0;
int * startR = rhsFact->R_etas_index+n3;
nCopyEnd=static_cast<int>((rhsFact->xeradr+nnetas)-startR);
nCopyStart=rhsFact->nnentu;
nCopyEnd = CoinMin(nCopyEnd+20,nnetas);
kCopyEnd = nnetas-nCopyEnd;
nCopyStart = CoinMin(nCopyStart+20,nnetas);
if (!n2&&!rhsFact->nnentu&&!rhsFact->nnentl) {
nCopyStart=nCopyEnd=0;
}
}
/* copy */
if(nCopyStart||nCopyEnd||true) {
#if 1
memcpy(fact->kw1adr,rhsFact->kw1adr,length*sizeof(double));
#else
c_ekkscpy((length*sizeof(double))/sizeof(int),
reinterpret_cast<int *>( rhsFact->kw1adr,reinterpret_cast<int *>( fact->kw1adr));
#endif
}
/* if malloc fails - we have lost memory - start again */
if (!fact->ndenuc &&fact->if_sparse_update) {
/* allow second copy of elements */
if (!canReuseEtas)
fact->xe2adr = clp_double(nnetas);
if (!fact->xe2adr) {
fact->maxNNetas=nnetas; /* dont allow any increase */
#ifdef PRINT_DEBUG
if (fact->if_sparse_update) {
printf("*** Sparse update off due to memory\n");
}
#endif
fact->if_sparse_update=0;
} else {
#ifndef NDEBUG
memset(fact->xe2adr,CLP_FILL,nnetas*sizeof(double));
#endif
}
} else {
clp_free(fact->xe2adr);
fact->xe2adr = 0;
fact->if_sparse_update=0;
}
if (!canReuseEtas)
fact->xeradr= clp_int(nnetas);
if (!fact->xeradr) {
nnetas=0;
} else {
#ifndef NDEBUG
memset(fact->xeradr,CLP_FILL,nnetas*sizeof(int));
#endif
/* copy */
if(nCopyStart||nCopyEnd) {
#if 0
memcpy(fact->xeradr,rhsFact->xeradr,nCopyStart*sizeof(int));
memcpy(fact->xeradr+kCopyEnd,rhsFact->xeradr+kCopyEnd,nCopyEnd*sizeof(int));
#else
c_ekkscpy(nCopyStart,rhsFact->xeradr,fact->xeradr);
c_ekkscpy(nCopyEnd,rhsFact->xeradr+kCopyEnd,fact->xeradr+kCopyEnd);
#endif
}
}
if (nnetas) {
if (!canReuseEtas)
fact->xecadr= clp_int(nnetas);
if (!fact->xecadr) {
nnetas=0;
} else {
#ifndef NDEBUG
memset(fact->xecadr,CLP_FILL,nnetas*sizeof(int));
#endif
/* copy */
if (fact->rows_ok&&(nCopyStart||nCopyEnd)) {
int i;
int * hcoliR=rhsFact->xecadr-1;
int * hcoli=fact->xecadr-1;
int * mrstrt=fact->xrsadr;
int * hinrow=fact->xrnadr;
#if 0
memcpy(fact->xecadr+kCopyEnd,rhsFact->xecadr+kCopyEnd,
nCopyEnd*sizeof(int));
#else
c_ekkscpy(nCopyEnd,rhsFact->xecadr+kCopyEnd,fact->xecadr+kCopyEnd);
#endif
if (!fact->xe2adr) {
for (i=0;i<fact->nrow;i++) {
int istart = mrstrt[i];
assert (istart>0&&istart<=nnetas);
assert (hinrow[i]>=0&&hinrow[i]<=fact->nrow);
memcpy(hcoli+istart,hcoliR+istart,hinrow[i]*sizeof(int));
}
} else {
double * de2valR=rhsFact->xe2adr-1;
double * de2val=fact->xe2adr-1;
#if 0
memcpy(fact->xe2adr+kCopyEnd,rhsFact->xe2adr+kCopyEnd,
nCopyEnd*sizeof(double));
#else
c_ekkdcpy(nCopyEnd,rhsFact->xe2adr+kCopyEnd
,fact->xe2adr+kCopyEnd);
#endif
for (i=0;i<fact->nrow;i++) {
int istart = mrstrt[i];
assert (istart>0&&istart<=nnetas);
assert (hinrow[i]>=0&&hinrow[i]<=fact->nrow);
memcpy(hcoli+istart,hcoliR+istart,hinrow[i]*sizeof(int));
memcpy(de2val+istart,de2valR+istart,hinrow[i]*sizeof(double));
#ifndef NDEBUG
{
int j;
for (j=istart;j<istart+hinrow[i];j++)
assert (fabs(de2val[j])<1.0e50);
}
#endif
}
}
}
}
}
if (nnetas) {
if (!canReuseEtas)
fact->xeeadr= clp_double(nnetas);
if (!fact->xeeadr) {
nnetas=0;
} else {
#ifndef NDEBUG
memset(fact->xeeadr,CLP_FILL,nnetas*sizeof(double));
#endif
/* copy */
if(nCopyStart||nCopyEnd) {
#if 0
memcpy(fact->xeeadr,rhsFact->xeeadr,nCopyStart*sizeof(double));
memcpy(fact->xeeadr+kCopyEnd,rhsFact->xeeadr+kCopyEnd,nCopyEnd*sizeof(double));
#else
c_ekkdcpy(nCopyStart,
rhsFact->xeeadr,fact->xeeadr);
c_ekkdcpy(nCopyEnd,
rhsFact->xeeadr+kCopyEnd,
fact->xeeadr+kCopyEnd);
#endif
}
/*fact->R_etas_index = &XERADR1()[kdnspt - 1];
fact->R_etas_element = &XEEADR1()[kdnspt - 1];*/
fact->R_etas_start = fact->xcsadr+
(rhsFact->R_etas_start-rhsFact->xcsadr);
fact->R_etas_index = fact->xeradr+
(rhsFact->R_etas_index-rhsFact->xeradr);
fact->R_etas_element = fact->xeeadr+
(rhsFact->R_etas_element-rhsFact->xeeadr);
}
}
}
assert (nnetas||!nrowmx);
fact->nnetas=nnetas;
clp_adjust_pointers(fact, -1);
clp_setup_pointers(fact);
clp_adjust_pointers(const_cast<EKKfactinfo *>(rhsFact), -1);
}
static void c_ekksmem_delete(EKKfactinfo *fact)
{
clp_adjust_pointers(fact, +1);
clp_free(fact->trueStart);
clp_free(fact->xe2adr);
clp_free(fact->xecadr);
clp_free(fact->xeradr);
clp_free(fact->xeeadr);
fact->eta_size=0;
fact->xrsadr = 0;
fact->xcsadr = 0;
fact->xrnadr = 0;
fact->xcnadr = 0;
fact->krpadr = 0;
fact->kcpadr = 0;
fact->xeradr = 0;
fact->xecadr = 0;
fact->xeeadr = 0;
fact->xe2adr = 0;
fact->trueStart = 0;
fact->kw2adr = 0;
fact->kw3adr = 0;
fact->kp1adr = 0;
fact->kp2adr = 0;
fact->kadrpm = 0;
fact->kw1adr = 0;
}
int c_ekk_IsSet(const int * array,int bit);
void c_ekk_Set(int * array,int bit);
void c_ekk_Unset(int * array,int bit);
int c_ekk_IsSet(const int * array,int bit)
{
int iWord = bit>>5;
int iBit = bit&31;
int word = array[iWord];
return (word&(1<<iBit))!=0;
}
void c_ekk_Set(int * array,int bit)
{
int iWord = bit>>5;
int iBit = bit&31;
array[iWord] |= (1<<iBit);
}
void c_ekk_Unset(int * array,int bit)
{
int iWord = bit>>5;
int iBit = bit&31;
array[iWord] &= ~(1<<iBit);
}
int CoinOslFactorization::factorize (
const CoinPackedMatrix & matrix,
int rowIsBasic[],
int columnIsBasic[],
double areaFactor )
{
setSolveMode(10);
if (areaFactor)
factInfo_.areaFactor = areaFactor;
const int * row = matrix.getIndices();
const CoinBigIndex * columnStart = matrix.getVectorStarts();
const int * columnLength = matrix.getVectorLengths();
const double * element = matrix.getElements();
int numberRows=matrix.getNumRows();
int numberColumns=matrix.getNumCols();
int numberBasic = 0;
CoinBigIndex numberElements=0;
int numberRowBasic=0;
// compute how much in basis
int i;
// Move pivot variables across if they look good
int * pivotTemp = new int [numberRows];
for (i=0;i<numberRows;i++) {
if (rowIsBasic[i]>=0)
pivotTemp[numberRowBasic++]=i;
}
numberBasic = numberRowBasic;
for (i=0;i<numberColumns;i++) {
if (columnIsBasic[i]>=0) {
pivotTemp[numberBasic++]=i;
numberElements += columnLength[i];
}
}
if ( numberBasic > numberRows ) {
return -2; // say too many in basis
}
numberElements = 3 * numberRows + 3 * numberElements + 20000;
setUsefulInformation(&numberRows,0);
getAreas ( numberRows, numberRows, numberElements,
2 * numberElements );
//fill
numberBasic=0;
numberElements=0;
// Fill in counts so we can skip part of preProcess
double * elementU=elements();
int * indexRowU=indices();
int * startColumnU=starts();
int * numberInRow=this->numberInRow();
int * numberInColumn=this->numberInColumn();
CoinZeroN ( numberInRow, numberRows );
CoinZeroN ( numberInColumn, numberRows );
for (i=0;i<numberRowBasic;i++) {
int iRow = pivotTemp[i];
// Change pivotTemp to correct sequence
pivotTemp[i]=iRow+numberColumns;
indexRowU[i]=iRow;
startColumnU[i]=i;
elementU[i]=-1.0;
numberInRow[iRow]=1;
numberInColumn[i]=1;
}
startColumnU[numberRowBasic]=numberRowBasic;
numberElements=numberRowBasic;
numberBasic=numberRowBasic;
for (i=0;i<numberColumns;i++) {
if (columnIsBasic[i]>=0) {
CoinBigIndex j;
for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
int iRow=row[j];
numberInRow[iRow]++;
indexRowU[numberElements]=iRow;
elementU[numberElements++]=element[j];
}
numberInColumn[numberBasic]=columnLength[i];
numberBasic++;
startColumnU[numberBasic]=numberElements;
}
}
preProcess ( );
factor ( );
if (status() == 0) {
int * pivotVariable = new int [numberRows];
postProcess(pivotTemp,pivotVariable);
for (i=0;i<numberRows;i++) {
int iPivot=pivotVariable[i];
if (iPivot<numberColumns) {
assert (columnIsBasic[iPivot]>=0);
columnIsBasic[iPivot]=i;
} else {
iPivot-=numberColumns;
assert (rowIsBasic[iPivot]>=0);
rowIsBasic[iPivot]=i;
}
}
delete [] pivotVariable;
}
delete [] pivotTemp;
return status_;
}
// Condition number - product of pivots after factorization
double
CoinOslFactorization::conditionNumber() const
{
double condition = 1.0;
const double *dluval = factInfo_.xeeadr+1-1; // stored before
const int *mcstrt = factInfo_.xcsadr+1;
for (int i=0;i<numberRows_;i++) {
const int kx = mcstrt[i];
const double dpiv = dluval[kx];
condition *= dpiv;
}
condition = CoinMax(fabs(condition),1.0e-50);
return 1.0/condition;
}