limp-cbc-0.3.2.0: cbits/coin/CoinOslFactorization3.cpp
/* $Id: CoinOslFactorization3.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 "CoinOslFactorization.hpp"
#include "CoinOslC.h"
#include "CoinFinite.hpp"
#define GO_DENSE 70
#define GO_DENSE_RATIO 1.8
int c_ekkclco(const EKKfactinfo *fact,int *hcoli,
int *mrstrt, int *hinrow, int xnewro);
void c_ekkclcp(const int *hcol, const double *dels, const int * mrstrt,
int *hrow, double *dels2, int *mcstrt,
int *hincol, int itype, int nnrow, int nncol,
int ninbas);
int c_ekkcmfc(EKKfactinfo *fact,
EKKHlink *rlink, EKKHlink *clink,
EKKHlink *mwork, void *maction_void,
int nnetas,
int *nsingp, int *xrejctp,
int *xnewrop, int xnewco,
int *ncompactionsp);
int c_ekkcmfy(EKKfactinfo *fact,
EKKHlink *rlink, EKKHlink *clink,
EKKHlink *mwork, void *maction_void,
int nnetas,
int *nsingp, int *xrejctp,
int *xnewrop, int xnewco,
int *ncompactionsp);
int c_ekkcmfd(EKKfactinfo *fact,
int *mcol,
EKKHlink *rlink, EKKHlink *clink,
int *maction,
int nnetas,
int *nnentlp, int *nnentup,
int *nsingp);
int c_ekkford(const EKKfactinfo *fact,const int *hinrow, const int *hincol,
int *hpivro, int *hpivco,
EKKHlink *rlink, EKKHlink *clink);
void c_ekkrowq(int *hrow, int *hcol, double *dels,
int *mrstrt,
const int *hinrow, int nnrow, int ninbas);
int c_ekkrwco(const EKKfactinfo *fact,double *dluval, int *hcoli, int *
mrstrt, int *hinrow, int xnewro);
int c_ekkrwcs(const EKKfactinfo *fact,double *dluval, int *hcoli, int *mrstrt,
const int *hinrow, const EKKHlink *mwork,
int nfirst);
void c_ekkrwct(const EKKfactinfo *fact,double *dluval, int *hcoli, int *mrstrt,
const int *hinrow, const EKKHlink *mwork,
const EKKHlink *rlink,
const short *msort, double *dsort,
int nlast, int xnewro);
int c_ekkshff(EKKfactinfo *fact,
EKKHlink *clink, EKKHlink *rlink,
int xnewro);
void c_ekkshfv(EKKfactinfo *fact, EKKHlink *rlink, EKKHlink *clink,
int xnewro);
int c_ekktria(EKKfactinfo *fact,
EKKHlink * rlink,
EKKHlink * clink,
int *nsingp,
int *xnewcop, int *xnewrop,
int *nlrowtp,
const int ninbas);
#if 0
static void c_ekkafpv(int *hentry, int *hcoli,
double *dluval, int *mrstrt,
int *hinrow, int nentry)
{
int j;
int nel, krs;
int koff;
int irow;
int ientry;
int * index;
for (ientry = 0; ientry < nentry; ++ientry) {
#ifdef INTEL
int * els_long,maxaij_long;
#endif
double * els;
irow = UNSHIFT_INDEX(hentry[ientry]);
nel = hinrow[irow];
krs = mrstrt[irow];
index=&hcoli[krs];
els=&dluval[krs];
#ifdef INTEL
els_long=reinterpret_cast<int *> (els);
maxaij_long=0;
#else
double maxaij = 0.f;
#endif
koff = 0;
j=0;
if ((nel&1)!=0) {
#ifdef INTEL
maxaij_long = els_long[1] & 0x7fffffff;
#else
maxaij=fabs(els[0]);
#endif
j=1;
}
while (j<nel) {
#ifdef INTEL
UNROLL_LOOP_BODY2({
int d_long = els_long[1+(j<<1)] & 0x7fffffff;
if (maxaij_long < d_long) {
maxaij_long = d_long;
koff=j;
}
j++;
});
#else
UNROLL_LOOP_BODY2({
double d = fabs(els[j]);
if (maxaij < d) {
maxaij = d;
koff=j;
}
j++;
});
#endif
}
SWAP(int, index[koff], index[0]);
SWAP(double, els[koff], els[0]);
}
} /* c_ekkafpv */
#endif
/* Uwe H. Suhl, March 1987 */
/* This routine processes col singletons during the LU-factorization. */
/* Return codes (checked version 1.11): */
/* 0: ok */
/* 6: pivot element too small */
/* -43: system error at label 420 (ipivot not found) */
/*
* This routine processes singleton columns during factorization of the
* nucleus. It is very similar to the first part of c_ekktria,
* but is more complex, because we now have to maintain the length
* lists.
* The differences are:
* (1) here we use the length list for length 1 rather than a queue.
* This routine is only called if it is known that there is a singleton
* column.
*
* (2) here we maintain hrowi by moving the last entry into the pivot
* column entry; that means we don't have to search for the pivot row
* entry like we do in c_ekktria.
*
* (3) here the hlink data structure is in use for the length lists,
* so we maintain it as we shorten rows and removing columns altogether.
*
*/
int c_ekkcsin(EKKfactinfo *fact,
EKKHlink *rlink, EKKHlink *clink,
int *nsingp)
{
#if 1
int *hcoli = fact->xecadr;
double *dluval = fact->xeeadr;
//double *dvalpv = fact->kw3adr;
int *mrstrt = fact->xrsadr;
int *hrowi = fact->xeradr;
int *mcstrt = fact->xcsadr;
int *hinrow = fact->xrnadr;
int *hincol = fact->xcnadr;
int *hpivro = fact->krpadr;
int *hpivco = fact->kcpadr;
#endif
const int nrow = fact->nrow;
const double drtpiv = fact->drtpiv;
int j, k, kc, kce, kcs, nzj;
double pivot;
int kipis, kipie;
int jpivot;
#ifndef NDEBUG
int kpivot=-1;
#else
int kpivot=-1;
#endif
bool small_pivot = false;
/* next singleton column.
* Note that when the pivot column itself was removed from the
* list, the column in the list after it (if any) moves to the
* head of the list.
* Also, if any column from the pivot row was reduced to length 1,
* then it will have been added to the list and now be in front.
*/
for (jpivot = hpivco[1]; jpivot > 0; jpivot = hpivco[1]) {
const int ipivot = hrowi[mcstrt[jpivot]]; /* (2) */
assert(ipivot);
/* The pivot row is being eliminated (3) */
C_EKK_REMOVE_LINK(hpivro, hinrow, rlink, ipivot);
/* Loop over nonzeros in pivot row: */
kipis = mrstrt[ipivot];
kipie = kipis + hinrow[ipivot] - 1;
for (k = kipis; k <= kipie; ++k) {
j = hcoli[k];
/*
* We're eliminating column jpivot,
* so we're eliminating the row it occurs in,
* so every column in this row is becoming one shorter.
*
* I don't know why we don't do the same for rejected columns.
*
* if xrejct is false, then no column has ever been rejected
* and this test wouldn't have to be made.
* However, that means this whole loop would have to be copied.
*/
if (! (clink[j].pre > nrow)) {
C_EKK_REMOVE_LINK(hpivco, hincol, clink, j); /* (3) */
}
--hincol[j];
kcs = mcstrt[j];
kce = kcs + hincol[j];
for (kc = kcs; kc <= kce; ++kc) {
if (ipivot == hrowi[kc]) {
break;
}
}
/* ASSERT !(kc>kce) */
/* (2) */
hrowi[kc] = hrowi[kce];
hrowi[kce] = 0;
if (j == jpivot) {
/* remember the slot corresponding to the pivot column */
kpivot = k;
}
else {
/*
* We just reduced the length of the column.
* If we haven't eliminated all of its elements completely,
* then we have to put it back in its new length list.
*
* If the column was rejected, we only put it back in a length
* list when it has been reduced to a singleton column,
* because it would just be rejected again.
*/
nzj = hincol[j];
if (! (nzj <= 0) &&
! (clink[j].pre > nrow && nzj != 1)) {
C_EKK_ADD_LINK(hpivco, nzj, clink, j); /* (3) */
}
}
}
assert (kpivot>0);
/* store pivot sequence number */
++fact->npivots;
rlink[ipivot].pre = -fact->npivots;
clink[jpivot].pre = -fact->npivots;
/* compute how much room we'll need later */
fact->nuspike += hinrow[ipivot];
/* check the pivot */
pivot = dluval[kpivot];
if (fabs(pivot) < drtpiv) {
/* pivot element too small */
small_pivot = true;
rlink[ipivot].pre = -nrow - 1;
clink[jpivot].pre = -nrow - 1;
++(*nsingp);
}
/* swap the pivoted column entry with the first entry in the row */
dluval[kpivot] = dluval[kipis];
dluval[kipis] = pivot;
hcoli[kpivot] = hcoli[kipis];
hcoli[kipis] = jpivot;
}
return (small_pivot);
} /* c_ekkcsin */
/* Uwe H. Suhl, March 1987 */
/* This routine processes row singletons during the computation of */
/* an LU-decomposition for the nucleus. */
/* Return codes (checked version 1.16): */
/* -5: not enough space in row file */
/* -6: not enough space in column file */
/* 7: pivot element too small */
/* -52: system error at label 220 (ipivot not found) */
/* -53: system error at label 400 (jpivot not found) */
int c_ekkrsin(EKKfactinfo *fact,
EKKHlink *rlink, EKKHlink *clink,
EKKHlink *mwork, int nfirst,
int *nsingp,
int *xnewcop, int *xnewrop,
int *nnentup,
int *kmxetap, int *ncompactionsp,
int *nnentlp)
{
#if 1
int *hcoli = fact->xecadr;
double *dluval = fact->xeeadr;
//double *dvalpv = fact->kw3adr;
int *mrstrt = fact->xrsadr;
int *hrowi = fact->xeradr;
int *mcstrt = fact->xcsadr;
int *hinrow = fact->xrnadr;
int *hincol = fact->xcnadr;
int *hpivro = fact->krpadr;
int *hpivco = fact->kcpadr;
#endif
const int nrow = fact->nrow;
const double drtpiv = fact->drtpiv;
int xnewro = *xnewrop;
int xnewco = *xnewcop;
int kmxeta = *kmxetap;
int nnentu = *nnentup;
int ncompactions = *ncompactionsp;
int nnentl = *nnentlp;
int i, j, k, kc, kr, npr, nzi;
double pivot;
int kjpis, kjpie, knprs, knpre;
double elemnt, maxaij;
int ipivot, epivco, lstart;
#ifndef NDEBUG
int kpivot=-1;
#else
int kpivot=-1;
#endif
int irtcod = 0;
const int nnetas = fact->nnetas;
lstart = nnetas - nnentl + 1;
for (ipivot = hpivro[1]; ipivot > 0; ipivot = hpivro[1]) {
const int jpivot = hcoli[mrstrt[ipivot]];
kjpis = mcstrt[jpivot];
kjpie = kjpis + hincol[jpivot] ;
for (k = kjpis; k < kjpie; ++k) {
i = hrowi[k];
/*
* We're eliminating row ipivot,
* so we're eliminating the column it occurs in,
* so every row in this column is becoming one shorter.
*
* No exception is made for rejected rows.
*/
C_EKK_REMOVE_LINK(hpivro, hinrow, rlink, i);
}
/* The pivot column is being eliminated */
/* I don't know why there is an exception for rejected columns */
if (! (clink[jpivot].pre > nrow)) {
C_EKK_REMOVE_LINK(hpivco, hincol, clink, jpivot);
}
epivco = hincol[jpivot] - 1;
kjpie = kjpis + epivco;
for (kc = kjpis; kc <= kjpie; ++kc) {
if (ipivot == hrowi[kc]) {
break;
}
}
/* ASSERT !(kc>kjpie) */
/* move the last column entry into this deleted one to keep */
/* the entries compact */
hrowi[kc] = hrowi[kjpie];
hrowi[kjpie] = 0;
/* store pivot sequence number */
++fact->npivots;
rlink[ipivot].pre = -fact->npivots;
clink[jpivot].pre = -fact->npivots;
/* Check if row or column files have to be compressed */
if (! (xnewro + epivco < lstart)) {
if (! (nnentu + epivco < lstart)) {
return (-5);
}
{
int iput = c_ekkrwcs(fact,dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
kmxeta += xnewro - iput ;
xnewro = iput - 1;
++ncompactions;
}
}
if (! (xnewco + epivco < lstart)) {
if (! (nnentu + epivco < lstart)) {
return (-5);
}
xnewco = c_ekkclco(fact,hrowi, mcstrt, hincol, xnewco);
++ncompactions;
}
/* This column has no more entries in it */
hincol[jpivot] = 0;
/* Perform numerical part of elimination. */
pivot = dluval[mrstrt[ipivot]];
if (fabs(pivot) < drtpiv) {
irtcod = 7;
rlink[ipivot].pre = -nrow - 1;
clink[jpivot].pre = -nrow - 1;
++(*nsingp);
}
/* If epivco is 0, then we can treat this like a singleton column (?)*/
if (! (epivco <= 0)) {
++fact->xnetal;
mcstrt[fact->xnetal] = lstart - 1;
hpivco[fact->xnetal] = ipivot;
/* Loop over nonzeros in pivot column. */
kjpis = mcstrt[jpivot];
kjpie = kjpis + epivco ;
nnentl+=epivco;
nnentu-=epivco;
for (kc = kjpis; kc < kjpie; ++kc) {
npr = hrowi[kc];
/* zero out the row entries as we go along */
hrowi[kc] = 0;
/* each row in the column is getting shorter */
--hinrow[npr];
/* find the entry in this row for the pivot column */
knprs = mrstrt[npr];
knpre = knprs + hinrow[npr];
for (kr = knprs; kr <= knpre; ++kr) {
if (jpivot == hcoli[kr])
break;
}
/* ASSERT !(kr>knpre) */
elemnt = dluval[kr];
/* move the last pivot column entry into this one */
/* to keep entries compact */
dluval[kr] = dluval[knpre];
hcoli[kr] = hcoli[knpre];
/*
* c_ekkmltf put the largest entries in front, and
* we want to maintain that property.
* There is only a problem if we just pivoted out the first
* entry, and there is more than one entry in the list.
*/
if (! (kr != knprs || hinrow[npr] <= 1)) {
maxaij = 0.f;
for (k = knprs; k <= knpre; ++k) {
if (! (fabs(dluval[k]) <= maxaij)) {
maxaij = fabs(dluval[k]);
kpivot = k;
}
}
assert (kpivot>0);
maxaij = dluval[kpivot];
dluval[kpivot] = dluval[knprs];
dluval[knprs] = maxaij;
j = hcoli[kpivot];
hcoli[kpivot] = hcoli[knprs];
hcoli[knprs] = j;
}
/* store elementary row transformation */
--lstart;
dluval[lstart] = -elemnt / pivot;
hrowi[lstart] = SHIFT_INDEX(npr);
/* Only add the row back in a length list if it isn't empty */
nzi = hinrow[npr];
if (! (nzi <= 0)) {
C_EKK_ADD_LINK(hpivro, nzi, rlink, npr);
}
}
++fact->nuspike;
}
}
*xnewrop = xnewro;
*xnewcop = xnewco;
*kmxetap = kmxeta;
*nnentup = nnentu;
*ncompactionsp = ncompactions;
*nnentlp = nnentl;
return (irtcod);
} /* c_ekkrsin */
int c_ekkfpvt(const EKKfactinfo *fact,
EKKHlink *rlink, EKKHlink *clink,
int *nsingp, int *xrejctp,
int *xipivtp, int *xjpivtp)
{
double zpivlu = fact->zpivlu;
#if 1
int *hcoli = fact->xecadr;
double *dluval = fact->xeeadr;
//double *dvalpv = fact->kw3adr;
int *mrstrt = fact->xrsadr;
int *hrowi = fact->xeradr;
int *mcstrt = fact->xcsadr;
int *hinrow = fact->xrnadr;
int *hincol = fact->xcnadr;
int *hpivro = fact->krpadr;
int *hpivco = fact->kcpadr;
#endif
int i, j, k, ke, kk, ks, nz, nz1, kce, kcs, kre, krs;
double minsze;
int marcst, mincst, mincnt, trials, nentri;
int jpivot=-1;
bool rjectd;
int ipivot;
const int nrow = fact->nrow;
int irtcod = 0;
/* this used to be initialized in c_ekklfct */
const int xtrial = 1;
trials = 0;
ipivot = 0;
mincst = COIN_INT_MAX;
mincnt = COIN_INT_MAX;
for (nz = 2; nz <= nrow; ++nz) {
nz1 = nz - 1;
if (mincnt <= nz) {
goto L900;
}
/* Search rows for a pivot */
for (i = hpivro[nz]; ! (i <= 0); i = rlink[i].suc) {
ks = mrstrt[i];
ke = ks + nz - 1;
/* Determine magnitude of minimal acceptable element */
minsze = fabs(dluval[ks]) * zpivlu;
for (k = ks; k <= ke; ++k) {
/* Consider a column only if it passes the stability test */
if (! (fabs(dluval[k]) < minsze)) {
j = hcoli[k];
marcst = nz1 * hincol[j];
if (! (marcst >= mincst)) {
mincst = marcst;
mincnt = hincol[j];
ipivot = i;
jpivot = j;
if (mincnt <= nz + 1) {
goto L900;
}
}
}
}
++trials;
if (trials >= xtrial) {
goto L900;
}
}
/* Search columns for a pivot */
j = hpivco[nz];
while (! (j <= 0)) {
/* XSEARD = XSEARD + 1 */
rjectd = false;
kcs = mcstrt[j];
kce = kcs + nz - 1;
for (k = kcs; k <= kce; ++k) {
i = hrowi[k];
nentri = hinrow[i];
marcst = nz1 * nentri;
if (! (marcst >= mincst)) {
/* Determine magnitude of minimal acceptable element */
minsze = fabs(dluval[mrstrt[i]]) * zpivlu;
krs = mrstrt[i];
kre = krs + nentri - 1;
for (kk = krs; kk <= kre; ++kk) {
if (hcoli[kk] == j)
break;
}
/* ASSERT (kk <= kre) */
/* perform stability test */
if (! (fabs(dluval[kk]) < minsze)) {
mincst = marcst;
mincnt = nentri;
ipivot = i;
jpivot = j;
rjectd = false;
if (mincnt <= nz) {
goto L900;
}
}
else {
if (ipivot == 0) {
rjectd = true;
}
}
}
}
++trials;
if (trials >= xtrial && ipivot > 0) {
goto L900;
}
if (rjectd) {
int jsuc = clink[j].suc;
++(*xrejctp);
C_EKK_REMOVE_LINK(hpivco, hincol, clink, j);
clink[j].pre = nrow + 1;
j = jsuc;
}
else {
j = clink[j].suc;
}
}
}
/* FLAG REJECTED ROWS (should this be columns ?) */
for (j = 1; j <= nrow; ++j) {
if (hinrow[j] == 0) {
rlink[j].pre = -nrow - 1;
++(*nsingp);
}
}
irtcod = 10;
L900:
*xipivtp = ipivot;
*xjpivtp = jpivot;
return (irtcod);
} /* c_ekkfpvt */
void c_ekkprpv(EKKfactinfo *fact,
EKKHlink *rlink, EKKHlink *clink,
int xrejct,
int ipivot, int jpivot)
{
#if 1
int *hcoli = fact->xecadr;
double *dluval = fact->xeeadr;
//double *dvalpv = fact->kw3adr;
int *mrstrt = fact->xrsadr;
int *hrowi = fact->xeradr;
int *mcstrt = fact->xcsadr;
int *hinrow = fact->xrnadr;
int *hincol = fact->xcnadr;
int *hpivro = fact->krpadr;
int *hpivco = fact->kcpadr;
#endif
int i, k;
int kc;
double pivot;
int kipis = mrstrt[ipivot];
int kipie = kipis + hinrow[ipivot] - 1;
#ifndef NDEBUG
int kpivot=-1;
#else
int kpivot=-1;
#endif
const int nrow = fact->nrow;
/* Update data structures */
{
int kjpis = mcstrt[jpivot];
int kjpie = kjpis + hincol[jpivot] ;
for (k = kjpis; k < kjpie; ++k) {
i = hrowi[k];
C_EKK_REMOVE_LINK(hpivro, hinrow, rlink, i);
}
}
for (k = kipis; k <= kipie; ++k) {
int j = hcoli[k];
if ((xrejct == 0) ||
! (clink[j].pre > nrow)) {
C_EKK_REMOVE_LINK(hpivco, hincol, clink, j);
}
--hincol[j];
int kcs = mcstrt[j];
int kce = kcs + hincol[j];
for (kc = kcs; kc < kce ; kc ++) {
if (hrowi[kc] == ipivot)
break;
}
assert (kc<kce||hrowi[kce]==ipivot);
hrowi[kc] = hrowi[kce];
hrowi[kce] = 0;
if (j == jpivot) {
kpivot = k;
}
}
assert (kpivot>0);
/* Store the pivot sequence number */
++fact->npivots;
rlink[ipivot].pre = -fact->npivots;
clink[jpivot].pre = -fact->npivots;
pivot = dluval[kpivot];
dluval[kpivot] = dluval[kipis];
dluval[kipis] = pivot;
hcoli[kpivot] = hcoli[kipis];
hcoli[kipis] = jpivot;
} /* c_ekkprpv */
/*
* c_ekkclco is almost exactly like c_ekkrwco.
*/
int c_ekkclco(const EKKfactinfo *fact,int *hcoli, int *mrstrt, int *hinrow, int xnewro)
{
#if 0
int *hcoli = fact->xecadr;
double *dluval = fact->xeeadr;
double *dvalpv = fact->kw3adr;
int *mrstrt = fact->xrsadr;
int *hrowi = fact->xeradr;
int *mcstrt = fact->xcsadr;
int *hinrow = fact->xrnadr;
int *hincol = fact->xcnadr;
int *hpivro = fact->krpadr;
int *hpivco = fact->kcpadr;
#endif
int i, k, nz, kold;
int kstart;
const int nrow = fact->nrow;
for (i = 1; i <= nrow; ++i) {
nz = hinrow[i];
if (0 < nz) {
/* save the last column entry of row i in hinrow */
/* and replace that entry with -i */
k = mrstrt[i] + nz - 1;
hinrow[i] = hcoli[k];
hcoli[k] = -i;
}
}
kstart = 0;
kold = 0;
for (k = 1; k <= xnewro; ++k) {
if (hcoli[k] != 0) {
++kstart;
/* if this is the last entry for the row... */
if (hcoli[k] < 0) {
/* restore the entry */
i = -hcoli[k];
hcoli[k] = hinrow[i];
/* update mrstart and hinrow */
mrstrt[i] = kold + 1;
hinrow[i] = kstart - kold;
kold = kstart;
}
hcoli[kstart] = hcoli[k];
}
}
/* INSERTED INCASE CALLED FROM YTRIAN JJHF */
mrstrt[nrow + 1] = kstart + 1;
return (kstart);
} /* c_ekkclco */
#undef MACTION_T
#define COIN_OSL_CMFC
#define MACTION_T short int
int c_ekkcmfc(EKKfactinfo *fact,
EKKHlink *rlink, EKKHlink *clink,
EKKHlink *mwork, void *maction_void,
int nnetas,
int *nsingp, int *xrejctp,
int *xnewrop, int xnewco,
int *ncompactionsp)
#include "CoinOslC.h"
#undef COIN_OSL_CMFC
#undef MACTION_T
static int c_ekkidmx(int n, const double *dx)
{
int ret_val;
int i;
double dmax;
--dx;
/* Function Body */
if (n < 1) {
return (0);
}
if (n == 1) {
return (1);
}
ret_val = 1;
dmax = fabs(dx[1]);
for (i = 2; i <= n; ++i) {
if (fabs(dx[i]) > dmax) {
ret_val = i;
dmax = fabs(dx[i]);
}
}
return ret_val;
} /* c_ekkidmx */
/* Return codes in IRTCOD/IRTCOD are */
/* 4: numerical problems */
/* 5: not enough space in row file */
/* 6: not enough space in column file */
int c_ekkcmfd(EKKfactinfo *fact,
int *mcol,
EKKHlink *rlink, EKKHlink *clink,
int *maction,
int nnetas,
int *nnentlp, int *nnentup,
int *nsingp)
{
int *hcoli = fact->xecadr;
double *dluval = fact->xeeadr;
int *mrstrt = fact->xrsadr;
int *hrowi = fact->xeradr;
int *mcstrt = fact->xcsadr;
int *hinrow = fact->xrnadr;
int *hincol = fact->xcnadr;
int *hpivro = fact->krpadr;
int *hpivco = fact->kcpadr;
int nnentl = *nnentlp;
int nnentu = *nnentup;
int storeZero = fact->ndenuc;
int mkrs[8];
double dpivyy[8];
/* Local variables */
int i, j;
double d0, dx;
int nz, ndo, krs;
int kend, jcol;
int irow, iput, jrow, krxs;
int mjcol[8];
double pivot;
int count;
int ilast, isort;
double dpivx, dsave;
double dpivxx[8];
double multip;
int lstart, ndense, krlast, kcount, idense, ipivot,
jdense, kchunk, jpivot;
const int nrow = fact->nrow;
int irtcod = 0;
lstart = nnetas - nnentl + 1;
/* put list of columns in last HROWI */
/* fix row order once for all */
ndense = nrow - fact->npivots;
iput = ndense + 1;
for (i = 1; i <= nrow; ++i) {
if (hpivro[i] > 0) {
irow = hpivro[i];
for (j = 1; j <= nrow; ++j) {
--iput;
maction[iput] = irow;
irow = rlink[irow].suc;
if (irow == 0) {
break;
}
}
}
}
if (iput != 1) {
++(*nsingp);
}
else {
/* Use HCOLI just for last row */
ilast = maction[1];
krlast = mrstrt[ilast];
/* put list of columns in last HCOLI */
iput = 0;
for (i = 1; i <= nrow; ++i) {
if (clink[i].pre >= 0) {
hcoli[krlast + iput] = i;
++iput;
}
}
if (iput != ndense) {
++(*nsingp);
}
else {
ndo = ndense / 8;
/* do most */
for (kcount = 1; kcount <= ndo; ++kcount) {
idense = ndense;
isort = 8;
for (count = ndense; count >= ndense - 7; --count) {
ipivot = maction[count];
krs = mrstrt[ipivot];
--isort;
mkrs[isort] = krs;
}
isort = 8;
for (count = ndense; count >= ndense - 7; --count) {
/* Find a pivot element */
--isort;
ipivot = maction[count];
krs = mkrs[isort];
jcol = c_ekkidmx(idense, &dluval[krs]) - 1;
pivot = dluval[krs + jcol];
--idense;
mcol[count] = jcol;
mjcol[isort] = mcol[count];
dluval[krs + jcol] = dluval[krs + idense];
if (fabs(pivot) < fact->zeroTolerance) {
pivot = 0.;
dpivx = 0.;
} else {
dpivx = 1. / pivot;
}
dluval[krs + idense] = pivot;
dpivxx[isort] = dpivx;
for (j = isort - 1; j >= 0; --j) {
krxs = mkrs[j];
multip = -dluval[krxs + jcol] * dpivx;
dluval[krxs + jcol] = dluval[krxs + idense];
/* for moment skip if zero */
if (fabs(multip) > fact->zeroTolerance) {
for (i = 0; i < idense; ++i) {
dluval[krxs + i] += multip * dluval[krs + i];
}
} else {
multip = 0.;
}
dluval[krxs + idense] = multip;
}
}
/* sort all U in rows already done */
for (i = 7; i >= 0; --i) {
/* **** this is important bit */
krs = mkrs[i];
for (j = i - 1; j >= 0; --j) {
jcol = mjcol[j];
dsave = dluval[krs + jcol];
dluval[krs + jcol] = dluval[krs + idense + j];
dluval[krs + idense + j] = dsave;
}
}
/* leave IDENSE as it is */
if (ndense <= 400) {
for (jrow = ndense - 8; jrow >= 1; --jrow) {
irow = maction[jrow];
krxs = mrstrt[irow];
for (j = 7; j >= 0; --j) {
jcol = mjcol[j];
dsave = dluval[krxs + jcol];
dluval[krxs + jcol] = dluval[krxs + idense + j];
dluval[krxs + idense + j] = dsave;
}
for (j = 7; j >= 0; --j) {
krs = mkrs[j];
jdense = idense + j;
dpivx = dpivxx[j];
multip = -dluval[krxs + jdense] * dpivx;
if (fabs(multip) <= fact->zeroTolerance) {
multip = 0.;
}
dpivyy[j] = multip;
dluval[krxs + jdense] = multip;
for (i = idense; i < jdense; ++i) {
dluval[krxs + i] += multip * dluval[krs + i];
}
}
for (i = 0; i < idense; ++i) {
dx = dluval[krxs + i];
d0 = dpivyy[0] * dluval[mkrs[0] + i];
dx += dpivyy[1] * dluval[mkrs[1] + i];
d0 += dpivyy[2] * dluval[mkrs[2] + i];
dx += dpivyy[3] * dluval[mkrs[3] + i];
d0 += dpivyy[4] * dluval[mkrs[4] + i];
dx += dpivyy[5] * dluval[mkrs[5] + i];
d0 += dpivyy[6] * dluval[mkrs[6] + i];
dx += dpivyy[7] * dluval[mkrs[7] + i];
dluval[krxs + i] = d0 + dx;
}
}
} else {
for (jrow = ndense - 8; jrow >= 1; --jrow) {
irow = maction[jrow];
krxs = mrstrt[irow];
for (j = 7; j >= 0; --j) {
jcol = mjcol[j];
dsave = dluval[krxs + jcol];
dluval[krxs + jcol] = dluval[krxs + idense + j];
dluval[krxs + idense + j] = dsave;
}
for (j = 7; j >= 0; --j) {
krs = mkrs[j];
jdense = idense + j;
dpivx = dpivxx[j];
multip = -dluval[krxs + jdense] * dpivx;
if (fabs(multip) <= fact->zeroTolerance) {
multip = 0.;
}
dluval[krxs + jdense] = multip;
for (i = idense; i < jdense; ++i) {
dluval[krxs + i] += multip * dluval[krs + i];
}
}
}
for (kchunk = 0; kchunk < idense; kchunk += 400) {
kend = CoinMin(idense - 1, kchunk + 399);
for (jrow = ndense - 8; jrow >= 1; --jrow) {
irow = maction[jrow];
krxs = mrstrt[irow];
for (j = 7; j >= 0; --j) {
dpivyy[j] = dluval[krxs + idense + j];
}
for (i = kchunk; i <= kend; ++i) {
dx = dluval[krxs + i];
d0 = dpivyy[0] * dluval[mkrs[0] + i];
dx += dpivyy[1] * dluval[mkrs[1] + i];
d0 += dpivyy[2] * dluval[mkrs[2] + i];
dx += dpivyy[3] * dluval[mkrs[3] + i];
d0 += dpivyy[4] * dluval[mkrs[4] + i];
dx += dpivyy[5] * dluval[mkrs[5] + i];
d0 += dpivyy[6] * dluval[mkrs[6] + i];
dx += dpivyy[7] * dluval[mkrs[7] + i];
dluval[krxs + i] = d0 + dx;
}
}
}
}
/* resort all U in rows already done */
for (i = 7; i >= 0; --i) {
krs = mkrs[i];
for (j = 0; j < i; ++j) {
jcol = mjcol[j];
dsave = dluval[krs + jcol];
dluval[krs + jcol] = dluval[krs + idense + j];
dluval[krs + idense + j] = dsave;
}
}
ndense += -8;
}
idense = ndense;
/* do remainder */
for (count = ndense; count >= 1; --count) {
/* Find a pivot element */
ipivot = maction[count];
krs = mrstrt[ipivot];
jcol = c_ekkidmx(idense, &dluval[krs]) - 1;
pivot = dluval[krs + jcol];
--idense;
mcol[count] = jcol;
dluval[krs + jcol] = dluval[krs + idense];
if (fabs(pivot) < fact->zeroTolerance) {
dluval[krs + idense] = 0.;
} else {
dpivx = 1. / pivot;
dluval[krs + idense] = pivot;
for (jrow = idense; jrow >= 1; --jrow) {
irow = maction[jrow];
krxs = mrstrt[irow];
multip = -dluval[krxs + jcol] * dpivx;
dluval[krxs + jcol] = dluval[krxs + idense];
/* for moment skip if zero */
if (fabs(multip) > fact->zeroTolerance) {
dluval[krxs + idense] = multip;
for (i = 0; i < idense; ++i) {
dluval[krxs + i] += multip * dluval[krs + i];
}
} else {
dluval[krxs + idense] = 0.;
}
}
}
}
/* now create in form for OSL */
ndense = nrow - fact->npivots;
idense = ndense;
for (count = ndense; count >= 1; --count) {
/* Find a pivot element */
ipivot = maction[count];
krs = mrstrt[ipivot];
--idense;
jcol = mcol[count];
jpivot = hcoli[krlast + jcol];
++fact->npivots;
pivot = dluval[krs + idense];
if (pivot == 0.) {
hinrow[ipivot] = 0;
rlink[ipivot].pre = -nrow - 1;
++(*nsingp);
irtcod = 10;
} else {
rlink[ipivot].pre = -fact->npivots;
clink[jpivot].pre = -fact->npivots;
hincol[jpivot] = 0;
++fact->xnetal;
mcstrt[fact->xnetal] = lstart - 1;
hpivco[fact->xnetal] = ipivot;
for (jrow = idense; jrow >= 1; --jrow) {
irow = maction[jrow];
krxs = mrstrt[irow];
multip = dluval[krxs + idense];
/* for moment skip if zero */
if (multip != 0.||storeZero) {
/* Store elementary row transformation */
++nnentl;
--nnentu;
--lstart;
dluval[lstart] = multip;
hrowi[lstart] = SHIFT_INDEX(irow);
}
}
hcoli[krlast + jcol] = hcoli[krlast + idense];
/* update pivot row and last row HCOLI */
dluval[krs + idense] = dluval[krs];
hcoli[krlast + idense] = hcoli[krlast];
nz = 1;
dluval[krs] = pivot;
hcoli[krs] = jpivot;
if (!storeZero) {
for (i = 1; i <= idense; ++i) {
if (fabs(dluval[krs + i]) > fact->zeroTolerance) {
++nz;
hcoli[krs + nz - 1] = hcoli[krlast + i];
dluval[krs + nz - 1] = dluval[krs + i];
}
}
hinrow[ipivot] = nz;
} else {
for (i = 1; i <= idense; ++i) {
++nz;
hcoli[krs + nz - 1] = hcoli[krlast + i];
dluval[krs + nz - 1] = dluval[krs + i];
}
hinrow[ipivot] = nz;
}
}
}
}
}
*nnentlp = nnentl;
*nnentup = nnentu;
return (irtcod);
} /* c_ekkcmfd */
/* ***C_EKKCMFC */
/*
* Generate a variant of c_ekkcmfc that uses an maction array of type
* int rather than short.
*/
#undef MACTION_T
#define C_EKKCMFY
#define COIN_OSL_CMFC
#define MACTION_T int
int c_ekkcmfy(EKKfactinfo *fact,
EKKHlink *rlink, EKKHlink *clink,
EKKHlink *mwork, void *maction_void,
int nnetas,
int *nsingp, int *xrejctp,
int *xnewrop, int xnewco,
int *ncompactionsp)
#include "CoinOslC.h"
#undef COIN_OSL_CMFC
#undef C_EKKCMFY
#undef MACTION_T
int c_ekkford(const EKKfactinfo *fact,const int *hinrow, const int *hincol,
int *hpivro, int *hpivco,
EKKHlink *rlink, EKKHlink *clink)
{
int i, iri, nzi;
const int nrow = fact->nrow;
int nsing = 0;
/* Uwe H. Suhl, August 1986 */
/* Builds linked lists of rows and cols of nucleus for efficient */
/* pivot searching. */
memset(hpivro+1,0,nrow*sizeof(int));
memset(hpivco+1,0,nrow*sizeof(int));
for (i = 1; i <= nrow; ++i) {
//hpivro[i] = 0;
//hpivco[i] = 0;
assert(rlink[i].suc == 0);
assert(clink[i].suc == 0);
}
/* Generate double linked list of rows having equal numbers of */
/* nonzeros in each row. Skip pivotal rows. */
for (i = 1; i <= nrow; ++i) {
if (! (rlink[i].pre < 0)) {
nzi = hinrow[i];
if (nzi <= 0) {
++nsing;
rlink[i].pre = -nrow - 1;
}
else {
iri = hpivro[nzi];
hpivro[nzi] = i;
rlink[i].suc = iri;
rlink[i].pre = 0;
if (iri != 0) {
rlink[iri].pre = i;
}
}
}
}
/* Generate double linked list of cols having equal numbers of */
/* nonzeros in each col. Skip pivotal cols. */
for (i = 1; i <= nrow; ++i) {
if (! (clink[i].pre < 0)) {
nzi = hincol[i];
if (nzi <= 0) {
++nsing;
clink[i].pre = -nrow - 1;
}
else {
iri = hpivco[nzi];
hpivco[nzi] = i;
clink[i].suc = iri;
clink[i].pre = 0;
if (iri != 0) {
clink[iri].pre = i;
}
}
}
}
return (nsing);
} /* c_ekkford */
/* c version of OSL from 36100 */
/* Assumes that a basis exists in correct form */
/* Calls Uwe's routines (approximately) */
/* Then if OK shuffles U into column order */
/* Return codes: */
/* 0: everything ok */
/* 1: everything ok but performance would be better if more space */
/* would be make available */
/* 4: growth rate of element in U too big */
/* 5: not enough space in row file */
/* 6: not enough space in column file */
/* 7: pivot too small - col sing */
/* 8: pivot too small - row sing */
/* 10: matrix is singular */
/* I suspect c_ekklfct never returns 1 */
/*
* layout of data
*
* dluval/hcoli: (L^-1)B - hole - L factors
*
* The L factors are written from high to low, starting from nnetas.
* There are nnentl factors in L. lstart the next entry to use for the
* L factors. Eventually, (L^-1)B turns into U.
* The ninbas coefficients of matrix B are originally in the start of
* dluval/hcoli. As L transforms it, rows may have to be expanded.
* If there is room, they are copied to the start of the hole,
* otherwise the first part of this area is compacted, and hopefully
* there is then room.
* There are nnentu coefficients in (L^-1)B.
* nnentu + nnentl >= ninbas.
* nnentu + nnentl == ninbas if there has been no fill-in.
* nnentu is decreased when the pivot eliminates elements
* (in which case there is a corresponding increase in nnentl),
* and if pivoting happens to cancel out factors (in which case
* there is no corresponding increase in L).
* nnentu is increased if there is fill-in (no decrease in L).
* If nnentu + nnentl >= nnetas, then we've run out of room.
* It is not the case that the elements of (L^-1)B are all in the
* first nnentu positions of dluval/hcoli, but that is of course
* the lower bound on the number of positions needed to store it.
* nuspik is roughly the sum of the row lengths of the rows that were pivoted
* out. singleton rows in c_ekktria do not change nuspik, but
* c_ekkrsin does increment it for each singleton row.
* That is, there are nuspik elements that in the upper part of (L^-1)B,
* and (nnentu - nuspik) elements left in B.
*/
/*
* As part of factorization, we test candidate pivots for numerical
* stability; if the largest element in a row/col is much larger than
* the smallest, this generally causes problems. To easily determine
* what the largest element is, we ensure that it is always in front.
* This establishes this property; later on we take steps to preserve it.
*/
static void c_ekkmltf(const EKKfactinfo *fact,double *dluval, int *hcoli,
const int *mrstrt, const int *hinrow,
const EKKHlink *rlink)
{
#if 0
int *hcoli = fact->xecadr;
double *dluval = fact->xeeadr;
double *dvalpv = fact->kw3adr;
int *mrstrt = fact->xrsadr;
int *hrowi = fact->xeradr;
int *mcstrt = fact->xcsadr;
int *hinrow = fact->xrnadr;
int *hincol = fact->xcnadr;
int *hpivro = fact->krpadr;
int *hpivco = fact->kcpadr;
#endif
int i, j, k;
int koff=-1;
const int nrow = fact->nrow;
for (i = 1; i <= nrow; ++i) {
/* ignore rows that have already been pivoted */
/* if it is a singleton row, the property trivially holds */
if (! (rlink[i].pre < 0 || hinrow[i] <= 1)) {
const int krs = mrstrt[i];
const int kre = krs + hinrow[i] - 1;
double maxaij = 0.f;
/* this assumes that at least one of the dluvals is non-zero. */
for (k = krs; k <= kre; ++k) {
if (! (fabs(dluval[k]) <= maxaij)) {
maxaij = fabs(dluval[k]);
koff = k;
}
}
assert (koff>0);
maxaij = dluval[koff];
j = hcoli[koff];
dluval[koff] = dluval[krs];
hcoli[koff] = hcoli[krs];
dluval[krs] = maxaij;
hcoli[krs] = j;
}
}
} /* c_ekkmltf */
int c_ekklfct( register EKKfactinfo *fact)
{
const int nrow = fact->nrow;
int ninbas = fact->xcsadr[nrow+1]-1;
int ifvsol = fact->ifvsol;
int *hcoli = fact->xecadr;
double *dluval = fact->xeeadr;
int *mrstrt = fact->xrsadr;
int *hrowi = fact->xeradr;
int *mcstrt = fact->xcsadr;
int *hinrow = fact->xrnadr;
int *hincol = fact->xcnadr;
int *hpivro = fact->krpadr;
int *hpivco = fact->kcpadr;
EKKHlink *rlink = fact->kp1adr;
EKKHlink *clink = fact->kp2adr;
EKKHlink *mwork = (reinterpret_cast<EKKHlink*>(fact->kw1adr))-1;
int nsing, kdnspt, xnewro, xnewco;
int i;
int xrejct;
int irtcod;
const int nnetas = fact->nnetas;
int ncompactions;
double save_drtpiv = fact->drtpiv;
double save_zpivlu = fact->zpivlu;
if (ifvsol > 0 && fact->invok < 0) {
fact->zpivlu = CoinMin(0.9, fact->zpivlu * 10.);
fact->drtpiv=1.0e-8;
}
rlink --;
clink --;
/* Function Body */
hcoli[nnetas] = 1;
hrowi[nnetas] = 1;
dluval[nnetas] = 0.0;
/* set amount of work */
xrejct = 0;
nsing = 0;
kdnspt = nnetas + 1;
fact->ndenuc = 0;
/* Triangularize */
irtcod = c_ekktria(fact,rlink,clink,
&nsing,
&xnewco, &xnewro,
&ncompactions, ninbas);
fact->nnentl = ninbas - fact->nnentu;
if (irtcod < 0) {
/* no space or system error */
goto L8000;
}
if (irtcod != 0 && fact->invok >= 0) {
goto L8500; /* 7 or 8 - pivot too small */
}
#if 0
/* is this necessary ? */
lstart = nnetas - fact->nnentl + 1;
for (i = lstart; i <= nnetas; ++i) {
hrowi[i] = (hcoli[i] << 3);
}
#endif
/* See if finished */
if (! (fact->npivots >= nrow)) {
int nsing1;
/* No - do nucleus */
nsing1 = c_ekkford(fact,hinrow, hincol, hpivro, hpivco, rlink, clink);
nsing+= nsing1;
if (nsing1 != 0 && fact->invok >= 0) {
irtcod=7;
goto L8500;
}
c_ekkmltf(fact,dluval, hcoli, mrstrt, hinrow, rlink);
{
bool callcmfy = false;
if (nrow > 32767) {
int count = 0;
for (i = 1; i <= nrow; ++i) {
count = CoinMax(count,hinrow[i]);
}
if (count + nrow - fact->npivots > 32767) {
/* will have to use I*4 version of CMFC */
/* no changes to pointer params */
callcmfy = true;
}
}
irtcod = (callcmfy ? c_ekkcmfy : c_ekkcmfc)
(fact,
rlink, clink,
mwork, &mwork[nrow + 1],
nnetas,
&nsing, &xrejct,
&xnewro, xnewco,
&ncompactions);
/* irtcod one of 0,-5,7,10 */
}
if (irtcod < 0) {
goto L8000;
}
kdnspt = nnetas - fact->nnentl;
}
/* return if error */
if (nsing > 0 || irtcod == 10) {
irtcod = 99;
}
/* irtcod one of 0,7,99 */
if (irtcod != 0) {
goto L8500;
}
++fact->xnetal;
mcstrt[fact->xnetal] = nnetas - fact->nnentl;
/* give message if tight on memory */
if (ncompactions > 2 ) {
if (1) {
int etasize =CoinMax(4*fact->nnentu+(nnetas-fact->nnentl)+1000,fact->eta_size);
fact->eta_size=CoinMin(static_cast<int>(1.2*fact->eta_size),etasize);
if (fact->maxNNetas>0&&fact->eta_size>
fact->maxNNetas) {
fact->eta_size=fact->maxNNetas;
}
} /* endif */
}
/* Shuffle U and multiply L by 8 (if assembler) */
{
int jrtcod = c_ekkshff(fact, clink, rlink,
xnewro);
/* nR_etas is the number of R transforms;
* it is incremented only in c_ekketsj.
*/
fact->nR_etas = 0;
/*fact->R_etas_start = mcstrt+nrow+fact->nnentl+3;*/
fact->R_etas_start[1] = /*kdnspt - 1*/0; /* magic */
fact->R_etas_index = &fact->xeradr[kdnspt - 1];
fact->R_etas_element = &fact->xeeadr[kdnspt - 1];
if (jrtcod != 0) {
irtcod = jrtcod;
/* irtcod == 2 */
}
}
goto L8500;
/* Fatal error */
L8000:
if (1) {
if (fact->maxNNetas != fact->eta_size &&
nnetas) {
/* return and get more space */
/* double eta_size, unless that exceeds max (if there is one) */
fact->eta_size = fact->eta_size<<1;
if (fact->maxNNetas > 0 &&
fact->eta_size > fact->maxNNetas) {
fact->eta_size = fact->maxNNetas;
}
return (5);
}
}
/*c_ekkmesg_no_i1(121, -irtcod);*/
irtcod = 3;
L8500:
/* restore pivot tolerance */
fact->drtpiv=save_drtpiv;
fact->zpivlu=save_zpivlu;
#ifndef NDEBUG
if (fact->rows_ok) {
int * hinrow=fact->xrnadr;
if (!fact->xe2adr) {
for (int i=1;i<=fact->nrow;i++) {
assert (hinrow[i]>=0&&hinrow[i]<=fact->nrow);
}
}
}
#endif
return (irtcod);
} /* c_ekklfct */
/*
summary of return codes
c_ekktria:
7 small pivot
-5 no memory
c_ekkcsin:
returns true if small pivot
c_ekkrsin:
-5 no memory
7 small pivot
c_ekkfpvt:
10: no pivots found (singular)
c_ekkcmfd:
10: zero pivot (not just small)
c_ekkcmfc:
-5: no memory
any non-zero code from c_ekkcsin, c_ekkrsin, c_ekkfpvt, c_ekkprpv, c_ekkcmfd
c_ekkshff:
2: singular
c_ekklfct:
any positive code from c_ekktria, c_ekkcmfc, c_ekkshff (2,7,10)
*except* 10, which is changed to 99.
all negative return codes are changed to 5 or 3
(5 == ran out of memory but could get more,
3 == ran out of memory, no luck)
so: 2,3,5,7,99
c_ekklfct1:
1: c_ekksmem_invert failed
2: c_ekkslcf/c_ekkslct ran out of room
any return code from c_ekklfct, except 2 and 5
*/
void c_ekkrowq(int *hrow, int *hcol, double *dels,
int *mrstrt,
const int *hinrow, int nnrow, int ninbas)
{
int i, k, iak, jak;
double daik;
int iloc;
double dsave;
int isave, jsave;
/* Order matrix rowwise using MRSTRT, DELS, HCOL */
k = 1;
/* POSITION AFTER END OF ROW */
for (i = 1; i <= nnrow; ++i) {
k += hinrow[i];
mrstrt[i] = k;
}
for (k = ninbas; k >= 1; --k) {
iak = hrow[k];
if (iak != 0) {
daik = dels[k];
jak = hcol[k];
hrow[k] = 0;
while (1) {
--mrstrt[iak];
iloc = mrstrt[iak];
dsave = dels[iloc];
isave = hrow[iloc];
jsave = hcol[iloc];
dels[iloc] = daik;
hrow[iloc] = 0;
hcol[iloc] = jak;
if (isave == 0)
break;
daik = dsave;
iak = isave;
jak = jsave;
}
}
}
} /* c_ekkrowq */
int c_ekkrwco(const EKKfactinfo *fact,double *dluval,
int *hcoli, int *mrstrt, int *hinrow, int xnewro)
{
int i, k, nz, kold;
int kstart;
const int nrow = fact->nrow;
for (i = 1; i <= nrow; ++i) {
nz = hinrow[i];
if (0 < nz) {
/* save the last column entry of row i in hinrow */
/* and replace that entry with -i */
k = mrstrt[i] + nz - 1;
hinrow[i] = hcoli[k];
hcoli[k] = -i;
}
}
kstart = 0;
kold = 0;
for (k = 1; k <= xnewro; ++k) {
if (hcoli[k] != 0) {
++kstart;
/* if this is the last entry for the row... */
if (hcoli[k] < 0) {
/* restore the entry */
i = -hcoli[k];
hcoli[k] = hinrow[i];
/* update mrstart and hinrow */
/* ACTUALLY, hinrow should already be accurate */
mrstrt[i] = kold + 1;
hinrow[i] = kstart - kold;
kold = kstart;
}
/* move the entry */
dluval[kstart] = dluval[k];
hcoli[kstart] = hcoli[k];
}
}
return (kstart);
} /* c_ekkrwco */
int c_ekkrwcs(const EKKfactinfo *fact,double *dluval, int *hcoli, int *mrstrt,
const int *hinrow, const EKKHlink *mwork,
int nfirst)
{
#if 0
int *hcoli = fact->xecadr;
double *dluval = fact->xeeadr;
double *dvalpv = fact->kw3adr;
int *mrstrt = fact->xrsadr;
int *hrowi = fact->xeradr;
int *mcstrt = fact->xcsadr;
int *hinrow = fact->xrnadr;
int *hincol = fact->xcnadr;
int *hpivro = fact->krpadr;
int *hpivco = fact->kcpadr;
#endif
int i, k, k1, k2, nz;
int irow, iput;
const int nrow = fact->nrow;
/* Compress row file */
iput = 1;
irow = nfirst;
for (i = 1; i <= nrow; ++i) {
nz = hinrow[irow];
k1 = mrstrt[irow];
if (k1 != iput) {
mrstrt[irow] = iput;
k2 = k1 + nz - 1;
for (k = k1; k <= k2; ++k) {
dluval[iput] = dluval[k];
hcoli[iput] = hcoli[k];
++iput;
}
} else {
iput += nz;
}
irow = mwork[irow].suc;
}
return (iput);
} /* c_ekkrwcs */
void c_ekkrwct(const EKKfactinfo *fact,double *dluval, int *hcoli, int *mrstrt,
const int *hinrow, const EKKHlink *mwork,
const EKKHlink *rlink,
const short *msort, double *dsort,
int nlast, int xnewro)
{
#if 0
int *hcoli = fact->xecadr;
double *dluval = fact->xeeadr;
double *dvalpv = fact->kw3adr;
int *mrstrt = fact->xrsadr;
int *hrowi = fact->xeradr;
int *mcstrt = fact->xcsadr;
int *hinrow = fact->xrnadr;
int *hincol = fact->xcnadr;
int *hpivro = fact->krpadr;
int *hpivco = fact->kcpadr;
#endif
int i, k, k1, nz, icol;
int kmax;
int irow, iput;
int ilook;
const int nrow = fact->nrow;
iput = xnewro;
irow = nlast;
kmax = nrow - fact->npivots;
for (i = 1; i <= nrow; ++i) {
nz = hinrow[irow];
k1 = mrstrt[irow] - 1;
if (rlink[irow].pre < 0) {
/* pivoted on already */
iput -= nz;
if (k1 != iput) {
mrstrt[irow] = iput + 1;
for (k = nz; k >= 1; --k) {
dluval[iput + k] = dluval[k1 + k];
hcoli[iput + k] = hcoli[k1 + k];
}
}
} else {
/* not pivoted - going dense */
iput -= kmax;
mrstrt[irow] = iput + 1;
c_ekkdzero( kmax, &dsort[1]);
for (k = 1; k <= nz; ++k) {
icol = hcoli[k1 + k];
ilook = msort[icol];
dsort[ilook] = dluval[k1 + k];
}
c_ekkdcpy(kmax,
(dsort+1), (dluval+iput + 1));
}
irow = mwork[irow].pre;
}
} /* c_ekkrwct */
/* takes Uwe's modern structures and puts them back 20 years */
int c_ekkshff(EKKfactinfo *fact,
EKKHlink *clink, EKKHlink *rlink,
int xnewro)
{
int *hpivro = fact->krpadr;
int i, j;
int nbas, icol;
int ipiv;
const int nrow = fact->nrow;
int nsing;
for (i = 1; i <= nrow; ++i) {
j = -rlink[i].pre;
rlink[i].pre = j;
if (j > 0 && j <= nrow) {
hpivro[j] = i;
}
j = -clink[i].pre;
clink[i].pre = j;
}
/* hpivro[j] is now (hopefully) the row that was pivoted on step j */
/* rlink[i].pre is the step in which row i was pivoted */
nbas = 0;
nsing = 0;
/* Decide if permutation wanted */
fact->first_dense=nrow-fact->ndenuc+1+1;
fact->last_dense=nrow;
/* rlink[].suc is dead at this point */
/*
* replace the the basis index
* with the pivot (or permuted) index generated by factorization.
* This eventually goes into mpermu.
*/
for (icol = 1; icol <= nrow; ++icol) {
int ibasis = icol;
ipiv = clink[ibasis].pre;
if (0 < ipiv && ipiv <= nrow) {
rlink[ibasis].suc = ipiv;
++nbas;
}
}
nsing = nrow - nbas;
if (nsing > 0) {
abort();
}
/* if we reach here, then rlink[1..nrow].suc == clink[1..nrow].pre */
/* switch off sparse update if any dense section */
{
const int notMuchRoom = (fact->nnentu + xnewro + 10 > fact->nnetas - fact->nnentl);
/* must be same as in c_ekkshfv */
if (fact->ndenuc || notMuchRoom||nrow<C_EKK_GO_SPARSE) {
#if PRINT_DEBUG
if (fact->if_sparse_update) {
printf("**** Switching off sparse update - dense - c_ekkshff\n");
}
#endif
fact->if_sparse_update=0;
}
}
/* hpivro[1..nrow] is not read by c_ekkshfv */
c_ekkshfv(fact,
rlink, clink,
xnewro);
return (0);
} /* c_ekkshff */
/* sorts on indices dragging elements with */
static void c_ekk_sort2(int * key , double * array2,int number)
{
int minsize=10;
int n = number;
int sp;
int *v = key;
int *m, t;
int * ls[32] , * rs[32];
int *l , *r , c;
double it;
int j;
/*check already sorted */
#ifndef LONG_MAX
#define LONG_MAX 0x7fffffff;
#endif
int last=-LONG_MAX;
for (j=0;j<number;j++) {
if (key[j]>=last) {
last=key[j];
} else {
break;
} /* endif */
} /* endfor */
if (j==number) {
return;
} /* endif */
sp = 0 ; ls[sp] = v ; rs[sp] = v + (n-1) ;
while( sp >= 0 )
{
if ( rs[sp] - ls[sp] > minsize )
{
l = ls[sp] ; r = rs[sp] ; m = l + (r-l)/2 ;
if ( *l > *m )
{
t = *l ; *l = *m ; *m = t ;
it = array2[l-v] ; array2[l-v] = array2[m-v] ; array2[m-v] = it ;
}
if ( *m > *r )
{
t = *m ; *m = *r ; *r = t ;
it = array2[m-v] ; array2[m-v] = array2[r-v] ; array2[r-v] = it ;
if ( *l > *m )
{
t = *l ; *l = *m ; *m = t ;
it = array2[l-v] ; array2[l-v] = array2[m-v] ; array2[m-v] = it ;
}
}
c = *m ;
while ( r - l > 1 )
{
while ( *(++l) < c ) ;
while ( *(--r) > c ) ;
t = *l ; *l = *r ; *r = t ;
it = array2[l-v] ; array2[l-v] = array2[r-v] ; array2[r-v] = it ;
}
l = r - 1 ;
if ( l < m )
{ ls[sp+1] = ls[sp] ;
rs[sp+1] = l ;
ls[sp ] = r ;
}
else
{ ls[sp+1] = r ;
rs[sp+1] = rs[sp] ;
rs[sp ] = l ;
}
sp++ ;
}
else sp-- ;
}
for ( l = v , m = v + (n-1) ; l < m ; l++ )
{ if ( *l > *(l+1) )
{
c = *(l+1) ;
it = array2[(l-v)+1] ;
for ( r = l ; r >= v && *r > c ; r-- )
{
*(r+1) = *r ;
array2[(r-v)+1] = array2[(r-v)] ;
}
*(r+1) = c ;
array2[(r-v)+1] = it ;
}
}
}
/* For each row compute reciprocal of pivot element and take out of */
/* Also use HLINK(1 to permute column numbers */
/* and HPIVRO to permute row numbers */
/* Sort into column order as was stored by row */
/* If Assembler then shift row numbers in L by 3 */
/* Put column numbers in U for L-U update */
/* and multiply U elements by - reciprocal of pivot element */
/* and set up backward pointers for pivot rows */
void c_ekkshfv(EKKfactinfo *fact,
EKKHlink *rlink, EKKHlink *clink,
int xnewro)
{
int *hcoli = fact->xecadr;
double *dluval = fact->xeeadr;
double *dvalpv = fact->kw3adr;
int *mrstrt = fact->xrsadr;
int *hrowi = fact->xeradr;
int *mcstrt = fact->xcsadr;
int *hinrow = fact->xrnadr;
int *hincol = fact->xcnadr;
int *hpivro = fact->krpadr;
int *hpivco = fact->kcpadr;
double *dpermu = fact->kadrpm;
double * de2val = fact->xe2adr ? fact->xe2adr-1: 0;
int nnentu = fact->nnentu;
int xnetal = fact->xnetal;
int numberSlacks; /* numberSlacks not read */
int i, j, k, kk, nel;
int nroom;
bool need_more_space;
int ndenuc=fact->ndenuc;
int if_sparse_update=fact->if_sparse_update;
int nnentl = fact->nnentl;
int nnetas = fact->nnetas;
int *ihlink = (reinterpret_cast<int*> (clink))+1; /* can't use rlink for simple loop below */
const int nrow = fact->nrow;
const int maxinv = fact->maxinv;
/* this is not just a temporary - c_ekkbtrn etc use this */
int *mpermu = (reinterpret_cast<int*> (dpermu+nrow))+1;
int * temp = ihlink+nrow;
int * temp2 = temp+nrow;
const int notMuchRoom = (nnentu + xnewro + 10 > nnetas - nnentl);
/* compress hlink and make simpler */
for (i = 1; i <= nrow; ++i) {
mpermu[i] = rlink[i].pre;
ihlink[i] = rlink[i].suc;
}
/* mpermu[i] == the step in which row i was pivoted */
/* ihlink[i] == the step in which col i was pivoted */
/* must be same as in c_ekkshff */
if (fact->ndenuc||notMuchRoom||nrow<C_EKK_GO_SPARSE) {
int ninbas;
/* CHANGE COLUMN NUMBERS AND FILL IN RECIPROCALS */
/* ALSO RECOMPUTE NUMBER IN COLUMN */
/* initialize with a fake pivot in each column */
c_ekkscpy_0_1(nrow, 1, &hincol[1]);
if (notMuchRoom) {
fact->eta_size=static_cast<int>(1.05*fact->eta_size);
/* eta_size can be no larger than maxNNetas */
if (fact->maxNNetas > 0 &&
fact->eta_size > fact->maxNNetas) {
fact->eta_size=fact->maxNNetas;
}
} /* endif */
/* For each row compute reciprocal of pivot element and take out of U */
/* Also use ihlink to permute column numbers */
/* the rows are not stored compactly or in order,
* so we have to find out where the last one is stored */
ninbas=0;
for (i = 1; i <= nrow; ++i) {
int jpiv=mpermu[i];
int nin=hinrow[i];
int krs = mrstrt[i];
int kre = krs + nin;
temp[jpiv]=krs;
temp2[jpiv]=nin;
ninbas = CoinMax(kre, ninbas);
/* c_ekktria etc ensure that the first row entry is the pivot */
dvalpv[jpiv] = 1. / dluval[krs];
hcoli[krs] = 0; /* probably needed for c_ekkrowq */
/* room for the pivot has already been allocated, so hincol ok */
for (kk = krs + 1; kk < kre; ++kk) {
int j = ihlink[hcoli[kk]];
hcoli[kk] = j; /* permute the col index */
hrowi[kk] = jpiv; /* permute the row index */
++hincol[j];
}
}
/* temp [mpermu[i]] == mrstrt[i] */
/* temp2[mpermu[i]] == hinrow[i] */
ninbas--; /* ???? */
c_ekkscpy(nrow, &temp[1], &mrstrt[1]);
c_ekkscpy(nrow, &temp2[1], &hinrow[1]);
/* now mrstrt, hinrow, hcoli and hrowi have been permuted */
/* Sort into column order as was stored by row */
/* There will be an empty entry in front of each each column,
* because we initialized hincol to 1s, and c_ekkrowq fills in
* entries from the back */
c_ekkrowq(hcoli, hrowi, dluval, mcstrt, hincol, nrow, ninbas);
/* The shuffle zeroed out column pointers */
/* Put them back for L-U update */
/* Also multiply U elements by - reciprocal of pivot element */
/* Also decrement mcstrt/hincol to give "real" sizes */
for (i = 1; i <= nrow; ++i) {
int kx = --mcstrt[i];
nel = --hincol[i];
hrowi[kx] = nel;
dluval[kx] = dvalpv[i];
#ifndef NO_SHIFT
for (int j=kx+1;j<=kx+nel;j++)
hrowi[j] = SHIFT_INDEX(hrowi[j]);
#endif
}
/* sort dense part */
for (i=nrow-ndenuc+1; i<=nrow; i++) {
int kx = mcstrt[i]+1; /* "real" entries start after pivot */
int nel = hincol[i];
c_ekk_sort2(&hrowi[kx],&dluval[kx],nel);
}
/* Recompute number in U */
nnentu = mcstrt[nrow] + hincol[nrow];
mcstrt[nrow + 4] = nnentu + 1; /* magic - AND DEAD */
/* as not much room switch off fast etas */
mrstrt[1] = 0; /* magic */
fact->rows_ok = false;
i = nrow + maxinv + 5; /* DEAD */
} else {
/* *************************************** */
/* enough memory to do a bit faster */
/* For each row compute reciprocal of pivot element and */
/* take out of U */
/* Also use HLINK(1 to permute column numbers */
int ninbas=0;
int ilast; /* last available entry */
int spareSpace;
double * dluval2;
/*int * hlink2 = ihlink+nrow;
int * mrstrt2 = hlink2+nrow;*/
/* mwork has order of row copy */
EKKHlink *mwork = (reinterpret_cast<EKKHlink*>(fact->kw1adr))-1;
fact->rows_ok = true;
if (if_sparse_update) {
ilast=nnetas-nnentl;
} else {
/* missing out nnentl stuff */
ilast=nnetas;
}
spareSpace=ilast-nnentu;
need_more_space=false;
/* save clean row copy if enough room */
nroom = (spareSpace) / nrow;
if (nrow<10000) {
if (nroom < 10) {
need_more_space=true;
}
} else {
if (nroom < 5&&!if_sparse_update) {
need_more_space=true;
}
}
if (nroom > CoinMin(50,maxinv)) {
need_more_space=false;
}
if (need_more_space) {
if (if_sparse_update) {
int i1=fact->eta_size+10*nrow;
fact->eta_size=static_cast<int>(1.2*fact->eta_size);
if (i1>fact->eta_size) {
fact->eta_size=i1;
}
} else {
fact->eta_size=static_cast<int>(1.05*fact->eta_size);
}
} else {
if (nroom<11) {
if (if_sparse_update) {
int i1=fact->eta_size+(11-nroom)*nrow;
fact->eta_size=static_cast<int>(1.2*fact->eta_size);
if (i1>fact->eta_size) {
fact->eta_size=i1;
}
}
}
}
if (fact->maxNNetas>0&&fact->eta_size>
fact->maxNNetas) {
fact->eta_size=fact->maxNNetas;
}
{
/* we can swap de2val and dluval to save copying */
int * eta_last=mpermu+nrow*2+3;
int * eta_next=eta_last+nrow+2;
int last=0;
eta_last[0]=-1;
if (nnentl) {
/* went into c_ekkcmfc - if not then in order */
int next;
/*next=mwork[((nrow+1)<<1)+1];*/
next=mwork[nrow+1].pre;
#ifdef DEBUG
j=mrstrt[next];
#endif
for (i = 1; i <= nrow; ++i) {
int iperm=mpermu[next];
eta_next[last]=iperm;
eta_last[iperm]=last;
temp[iperm] = mrstrt[next];
temp2[iperm] = hinrow[next];
#ifdef DEBUG
if (mrstrt[next]!=j) abort();
j=mrstrt[next]+hinrow[next];
#endif
/*next= mwork[(next<<1)+2];*/
next= mwork[next].suc;
last=iperm;
}
} else {
#ifdef DEBUG
j=0;
#endif
for (i = 1; i <= nrow; ++i) {
int iperm=mpermu[i];
eta_next[last]=iperm;
eta_last[iperm]=last;
temp[iperm] = mrstrt[i];
temp2[iperm] = hinrow[i];
last=iperm;
#ifdef DEBUG
if (mrstrt[i]<=j) abort();
if (i>1&&mrstrt[i]!=j+hinrow[i-1]) abort();
j=mrstrt[i];
#endif
}
}
eta_next[last]=nrow+1;
eta_last[nrow+1]=last;
eta_next[nrow+1]=nrow+2;
c_ekkscpy(nrow, &temp[1], &mrstrt[1]);
c_ekkscpy(nrow, &temp2[1], &hinrow[1]);
i=eta_last[nrow+1];
ninbas=mrstrt[i]+hinrow[i]-1;
#ifdef DEBUG
if (spareSpace<ninbas) {
abort();
}
#endif
c_ekkizero( nrow, &hincol[1]);
#ifdef DEBUG
for (i=nrow; i>0; i--) {
int krs = mrstrt[i];
int jpiv = hcoli[krs];
if (ihlink[jpiv]!=i) abort();
}
#endif
for (i = 1; i <= ninbas; ++i) {
k = hcoli[i];
k = ihlink[k];
#ifdef DEBUG
if (k<=0||k>nrow) abort();
#endif
hcoli[i]=k;
hincol[k]++;
}
#ifdef DEBUG
for (i=nrow; i>0; i--) {
int krs = mrstrt[i];
int jpiv = hcoli[krs];
if (jpiv!=i) abort();
if (krs>ninbas) abort();
}
#endif
/* Sort into column order as was stored by row */
k = 1;
/* Position */
for (kk = 1; kk <= nrow; ++kk) {
nel=hincol[kk];
mcstrt[kk] = k;
hrowi[k]=nel-1;
k += hincol[kk];
hincol[kk]=0;
}
if (de2val) {
dluval2=de2val;
} else {
dluval2=dluval+ninbas;
}
nnentu = k-1;
mcstrt[nrow + 4] = nnentu + 1;
/* create column copy */
for (i=nrow; i>0; i--) {
int krs = mrstrt[i];
int kre = krs + hinrow[i];
hinrow[i]--;
mrstrt[i]++;
{
int kx = mcstrt[i];
/*nel = hincol[i];
if (hrowi[kx]!=nel) abort();
hrowi[kx] = nel-1;*/
dluval2[kx] = 1.0 /dluval[krs];
/*hincol[i]=0;*/
for (kk = krs + 1; kk < kre; ++kk) {
int j = hcoli[kk];
int iput = hincol[j]+1;
hincol[j]=iput;
iput+= mcstrt[j];
hrowi[iput] = SHIFT_INDEX(i);
dluval2[iput] = dluval[kk];
}
}
}
if (de2val) {
double * a=dluval;
double * address;
/* move first down */
i=eta_next[0];
{
int krs=mrstrt[i];
nel=hinrow[i];
for (j=1;j<=nel;j++) {
hcoli[j]=hcoli[j+krs-1];
dluval[j]=dluval[j+krs-1];
}
}
mrstrt[i]=1;
/****** swap dluval and de2val !!!! ******/
/* should work even for dspace */
/* move L part across */
address=fact->xeeadr+1;
fact->xeeadr=fact->xe2adr-1;
fact->xe2adr=address;
if (nnentl) {
int n=xnetal-nrow-maxinv-5;
int j1,j2;
int * mcstrt2=mcstrt+nrow+maxinv+4;
j2 = mcstrt2[1];
j1 = mcstrt2[n+1]+1;
#if 0
memcpy(de2val+j1,dluval+j1,(j2-j1+1)*sizeof(double));
#else
c_ekkdcpy(j2-j1+1,
(dluval+j1),(de2val+j1));
#endif
}
dluval = de2val;
de2val = a;
} else {
/* copy down dluval */
#if 0
memcpy(&dluval[1],&dluval2[1],ninbas*sizeof(double));
#else
c_ekkdcpy(ninbas,
(dluval2+1),(dluval+1));
#endif
}
/* sort dense part */
for (i=nrow-ndenuc+1;i<=nrow;i++) {
int kx = mcstrt[i]+1;
int nel = hincol[i];
c_ekk_sort2(&hrowi[kx],&dluval[kx],nel);
}
}
mrstrt[nrow + 1] = ilast + 1;
}
/* Find first non slack */
for (i = 1; i <= nrow; ++i) {
int kcs = mcstrt[i];
if (hincol[i] != 0 || dluval[kcs] != SLACK_VALUE) {
break;
}
}
numberSlacks = i - 1;
{
/* set slacks to 1 */
int * array = fact->krpadr + ( fact->nrowmx+2);
int nSet = (numberSlacks)>>5;
int n2 = (fact->nrowmx+32)>>5;
int i;
memset(array,0xff,nSet*sizeof(int));
memset(array+nSet,0,(n2-nSet)*sizeof(int));
for (i=nSet<<5;i<=numberSlacks;i++)
c_ekk_Set(array,i);
c_ekk_Unset(array,fact->nrow+1); /* make sure off end not slack */
#ifndef NDEBUG
for (i=1;i<=numberSlacks;i++)
assert (c_ekk_IsSet(array,i));
for (;i<=fact->nrow;i++)
assert (!c_ekk_IsSet(array,i));
#endif
}
/* and set up backward pointers */
/* clean up HPIVCO for fancy assembler stuff */
/* xnetal was initialized to nrow + maxinv + 4 in c_ekktria, and grows */
c_ekkscpy_0_1(maxinv + 1, 1, &hpivco[nrow+4]); /* magic */
hpivco[xnetal] = 1;
/* shuffle down for gaps so can get rid of hpivco for L */
{
const int lstart = nrow + maxinv + 5;
int n=xnetal-lstart ; /* number of L entries */
int add,iel;
int * hpivco_L = &hpivco[lstart];
int * mcstrt_L = &mcstrt[lstart];
if (nnentl) {
/* elements of L were stored in descending order in dluval/hcoli */
int kle = mcstrt_L[0];
int kls = mcstrt_L[n]+1;
if(if_sparse_update) {
int i2,iel;
int * mrstrt2 = &mrstrt[nrow];
/* need row copy of L */
/* hpivro is spare for counts; just used as a temp buffer */
c_ekkizero( nrow, &hpivro[1]);
/* permute L indices; count L row lengths */
for (iel = kls; iel <= kle; ++iel) {
int jrow = mpermu[UNSHIFT_INDEX(hrowi[iel])];
hpivro[jrow]++;
hrowi[iel] = SHIFT_INDEX(jrow);
}
{
int ibase=nnetas-nnentl+1;
int firstDoRow=0;
for (i=1;i<=nrow;i++) {
mrstrt2[i]=ibase;
if (hpivro[i]&&!firstDoRow) {
firstDoRow=i;
}
ibase+=hpivro[i];
hpivro[i]=mrstrt2[i];
}
if (!firstDoRow) {
firstDoRow=nrow+1;
}
mrstrt2[i]=ibase;
fact->firstDoRow = firstDoRow;
}
i2=mcstrt_L[n];
for (i = n-1; i >= 0; --i) {
int i1 = mcstrt_L[i];
int ipiv=hpivco_L[i];
ipiv=mpermu[ipiv];
hpivco_L[i]=ipiv;
for (iel=i2 ; iel < i1; iel++) {
int irow = UNSHIFT_INDEX(hrowi[iel+1]);
int iput=hpivro[irow];
hpivro[irow]=iput+1;
hcoli[iput]=ipiv;
de2val[iput]=dluval[iel+1];
}
i2=i1;
}
} else {
/* just permute row numbers */
for (j = 0; j < n; ++j) {
hpivco_L[j] = mpermu[hpivco_L[j]];
}
for (iel = kls; iel <= kle; ++iel) {
int jrow = mpermu[UNSHIFT_INDEX(hrowi[iel])];
hrowi[iel] = SHIFT_INDEX(jrow);
}
}
add=hpivco_L[n-1]-hpivco_L[0]-n+1;
if (add) {
int i;
int last = hpivco_L[n-1];
int laststart = mcstrt_L[n];
int base=hpivco_L[0]-1;
/* adjust so numbers match */
mcstrt_L-=base;
hpivco_L-=base;
mcstrt_L[last]=laststart;
for (i=n-1;i>=0;i--) {
int ipiv=hpivco_L[i+base];
while (ipiv<last) {
mcstrt_L[last-1]=laststart;
hpivco_L[last-1]=last;
last--;
}
laststart=mcstrt_L[i+base];
mcstrt_L[last-1]=laststart;
hpivco_L[last-1]=last;
last--;
}
xnetal+=add;
}
}
//int lstart=fact->lstart;
//const int * COIN_RESTRICT hpivco = fact->kcpadr;
fact->firstLRow = hpivco[lstart];
}
fact->nnentu = nnentu;
fact->xnetal = xnetal;
/* now we have xnetal * we can set up pointers */
clp_setup_pointers(fact);
/* this is the array used in c_ekkbtrn; it is passed to c_ekkbtju as hpivco.
* this gets modified by F-T as we pivot columns in and out.
*/
{
/* do new hpivco */
int * hpivco_new = fact->kcpadr+1;
int * back = &fact->kcpadr[2*nrow+maxinv+4];
/* set zeroth to stop illegal read */
back[0]=1;
hpivco_new[nrow+1]=nrow+1; /* deliberate loop for dense tests */
hpivco_new[0]=1;
for (i=1;i<=nrow;i++) {
hpivco_new[i]=i+1;
back[i+1]=i;
}
back[1]=0;
fact->first_dense = CoinMax(fact->first_dense,4);
fact->numberSlacks=numberSlacks;
fact->lastSlack=numberSlacks;
fact->firstNonSlack=hpivco_new[numberSlacks];
}
/* also zero out permute region and nonzero */
c_ekkdzero( nrow, (dpermu+1));
if (if_sparse_update) {
char * nonzero = reinterpret_cast<char *> (&mpermu[nrow+1]); /* used in c_ekkbtrn */
/*c_ekkizero(nrow,(int *)nonzero);*/
c_ekkczero(nrow,nonzero);
/*memset(nonzero,0,nrow*sizeof(int));*/ /* for faster method */
}
for (i = 1; i <= nrow; ++i) {
hpivro[mpermu[i]] = i;
}
} /* c_ekkshfv */
static void c_ekkclcp1(const int *hcol, const int * mrstrt,
int *hrow, int *mcstrt,
int *hincol, int nnrow, int nncol,
int ninbas)
{
int i, j, kc, kr, kre, krs, icol;
int iput;
/* Create columnwise storage of row indices */
kc = 1;
for (j = 1; j <= nncol; ++j) {
mcstrt[j] = kc;
kc += hincol[j];
hincol[j] = 0;
}
mcstrt[nncol + 1] = ninbas + 1;
for (i = 1; i <= nnrow; ++i) {
krs = mrstrt[i];
kre = mrstrt[i + 1] - 1;
for (kr = krs; kr <= kre; ++kr) {
icol = hcol[kr];
iput = hincol[icol];
hincol[icol] = iput + 1;
iput += mcstrt[icol];
hrow[iput] = i;
}
}
} /* c_ekkclcp */
inline void c_ekkclcp2(const int *hcol, const double *dels, const int * mrstrt,
int *hrow, double *dels2, int *mcstrt,
int *hincol, int nnrow, int nncol,
int ninbas)
{
int i, j, kc, kr, kre, krs, icol;
int iput;
/* Create columnwise storage of row indices */
kc = 1;
for (j = 1; j <= nncol; ++j) {
mcstrt[j] = kc;
kc += hincol[j];
hincol[j] = 0;
}
mcstrt[nncol + 1] = ninbas + 1;
for (i = 1; i <= nnrow; ++i) {
krs = mrstrt[i];
kre = mrstrt[i + 1] - 1;
for (kr = krs; kr <= kre; ++kr) {
icol = hcol[kr];
iput = hincol[icol];
hincol[icol] = iput + 1;
iput += mcstrt[icol];
hrow[iput] = i;
dels2[iput] = dels[kr];
}
}
} /* c_ekkclcp */
int c_ekkslcf( register const EKKfactinfo *fact)
{
int * hrow = fact->xeradr;
int * hcol = fact->xecadr;
double * dels = fact->xeeadr;
int * hinrow = fact->xrnadr;
int * hincol = fact->xcnadr;
int * mrstrt = fact->xrsadr;
int * mcstrt = fact->xcsadr;
const int nrow = fact->nrow;
int ninbas;
/* space for etas */
const int nnetas = fact->nnetas;
ninbas=mcstrt[nrow+1]-1;
/* Now sort */
if (ninbas << 1 > nnetas) {
/* Put it in row order */
int i,k;
c_ekkrowq(hrow, hcol, dels, mrstrt, hinrow, nrow, ninbas);
k = 1;
for (i = 1; i <= nrow; ++i) {
mrstrt[i] = k;
k += hinrow[i];
}
mrstrt[nrow + 1] = k;
/* make a column copy without the extra values */
c_ekkclcp1(hcol, mrstrt, hrow, mcstrt, hincol, nrow, nrow, ninbas);
} else {
/* Move elements up memory */
c_ekkdcpy(ninbas,
(dels+1), (dels+ninbas + 1));
/* make a row copy with the extra values */
c_ekkclcp2(hrow, &dels[ninbas], mcstrt, hcol, dels, mrstrt, hinrow, nrow, nrow, ninbas);
}
return (ninbas);
} /* c_ekkslcf */
/* Uwe H. Suhl, September 1986 */
/* Removes lower and upper triangular factors from the matrix. */
/* Code for routine: 102 */
/* Return codes: */
/* 0: ok */
/* -5: not enough space in row file */
/* 7: pivot too small - col sing */
/*
* This selects singleton columns and rows for the LU factorization.
* Singleton columns require no
*
* (1) Note that columns are processed using a queue, not a stack;
* this produces better pivots.
*
* (2) At most nrows elements are ever entered into the queue.
*
* (3) When pivoting singleton columns, every column that is part of
* the pivot row is shortened by one, including the singleton column
* itself; the hincol entries are updated appropriately.
* Thus, pivoting on a singleton column may create other singleton columns
* (but not singleton rows).
* The dual property is true for rows.
*
* (4) Row entries (hrowi) are not changed when pivoting singleton columns.
* Singleton columns that are created as a result of pivoting the
* rows of other singleton columns will therefore have row entries
* corresponding to those pivoted rows. Since we need to find the
* row entry for the row being pivoted, we have to
* search its row entries for the one whose hlink entry indicates
* that it has not yet been pivoted.
*
* (9) As a result of pivoting columns, sections in hrowi corresponding to
* pivoted columns are no longer needed, and entries in sections
* for non-pivoted columns may have entries corresponding to pivoted rows.
* This is why hrowi needs to be compacted.
*
* (5) When the row_pre and col_pre fields of the hlink struct contain
* negative values, they indicate that the row has been pivoted, and
* the negative of that value is the pivot order.
* That is the only use for these fields in this routine.
*
* (6) This routine assumes that hlink is initialized to zeroes.
* Under this assumption, the following is an invariant in this routine:
*
* (clink[i].pre < 0) ==> (hincol[i]==0)
*
* The converse is not true; see (15).
*
* The dual is also true, but only while pivoting singletong rows,
* since we don't update hinrow while pivoting columns;
* THESE VALUES ARE USED LATER, BUT I DON'T UNDERSTAND HOW YET.
*
* (7) hpivco is used for two purposes. The low end is used to implement the
* queue when pivoting columns; the high end is used to hold eta-matrix
* entries.
*
* (8) As a result of pivoting columns, for all i:1<=i<=nrow, either
* hinrow[i] has not changed
* or
* hinrow[i] = 0
* This is another way of saying that pivoting singleton columns cannot
* create singleton rows.
* The dual holds for hincol after pivoting rows.
*
* (10) In constrast to (4), while pivoting rows we
* do not let the hcoli get out-of-date. That is because as part of
* the process of numerical pivoting we have to find the row entries
* for all the rows in the pivot column, so we may as well keep the
* entries up to date. This is done by moving the last column entry
* for each row into the entry that was used for the pivot column.
*
* (11) When pivoting a column, we must find the pivot row entry in
* its row table. Sometimes we search for other things at the same time.
* The same is true for pivoting columns. This search should never
* fail.
*
* (12) Information concerning the eta matrices is stored in the high
* ends of arrays that are also used to store information concerning
* the basis; these arrays are: hpivco, mcstrt, dluval and hcoli.
* Information is only stored in these arrays as a part of pivoting
* singleton rows, since the only thing that needs to be saved as
* a part of pivoting singleton columns is which rows and columns were chosen,
* and this is stored in hlink.
* Since they have to share the same array, the eta information grows
* downward instead of upward. Eventually, eta information may grow
* down to the top of the basis information. As pivoting proceeds,
* more and more of this information is no longer needed, so when this
* happens we can try compacting the arrays to see if we can recover
* enough space. lstart points at the bottom entry in the arrays,
* xnewro/xnewco at the top of the basis information, and each time we
* pivot a singleton row we know that we will need exactly as many new
* entries as there are rows in the pivot column, so we can easily
* determine if we need more room. The variable maxinv may be used
* to reserve extra room when inversion starts.
*
* (13) Eta information is stored in a fashion that is similar to how
* matrices are stored. There is one entry in hpivco and mcstrt for
* each eta (other than the initial ones for singleton columns and
* for singleton rows that turn out to be singleton columns),
* in the order they were chosen. hpivco records the pivot row,
* and mcstrt points at the first entry in the other two arrays
* for this row. dluval contains the actual eta values for the column,
* and hcoli the rows these values were in.
* These entries in mcstrt and hpivco grow upward; they start above
* the entries used to store basis information.
* (Actually, I don't see why they need to start maxinv+4 entries past the top).
*
* (14) c_ekkrwco assumes that invalidated hrowi/hcoli entries contain 0.
*
* (15) When pivoting singleton columns, it may possibly happen
* that a row with all singleton column entries is created.
* In this case, all of the columns will be enqueued, and pivoting
* on any of them eliminates the rest, without their being chosen
* as pivots. The dual holds for singleton rows.
* DOES THIS INDICATE A SINGULARITY?
*
* (15) There are some aspects of the implementation that I find odd.
* hrowi is not set to 0 for pivot rows while pivoting singleton columns,
* which would make sense to me. Things don't work if this isn't done,
* so the information is used somehow later on. Also, the information
* for the pivot column is shifted to the front of the pivot row
* when pivoting singleton columns; this is also necessary for reasons
* I don't understand.
*/
int c_ekktria(EKKfactinfo *fact,
EKKHlink * rlink,
EKKHlink * clink,
int *nsingp,
int *xnewcop, int *xnewrop,
int *ncompactionsp,
const int ninbas)
{
const int nrow = fact->nrow;
const int maxinv = fact->maxinv;
int *hcoli = fact->xecadr;
double *dluval = fact->xeeadr;
int *mrstrt = fact->xrsadr;
int *hrowi = fact->xeradr;
int *mcstrt = fact->xcsadr;
int *hinrow = fact->xrnadr;
int *hincol = fact->xcnadr;
int *stack = fact->krpadr; /* normally hpivro */
int *hpivco = fact->kcpadr;
const double drtpiv = fact->drtpiv;
CoinZeroN(reinterpret_cast<int *>(rlink+1),static_cast<int>(nrow*(sizeof(EKKHlink)/sizeof(int))));
CoinZeroN(reinterpret_cast<int *>(clink+1),static_cast<int>(nrow*(sizeof(EKKHlink)/sizeof(int))));
fact->npivots = 0;
/* Use NUSPIK to keep sum of deactivated row counts */
fact->nuspike = 0;
int xnetal = nrow + maxinv + 4;
int xnewro = mrstrt[nrow] + hinrow[nrow] - 1;
int xnewco = xnewro;
int kmxeta = ninbas;
int ncompactions = 0;
int i, j, k, kc, kce, kcs, npr;
double pivot;
int kipis, kipie, kjpis, kjpie, knprs, knpre;
int ipivot, jpivot, stackc, stackr;
#ifndef NDEBUG
int kpivot=-1;
#else
int kpivot=-1;
#endif
int epivco, kstart, maxstk;
int irtcod = 0;
int lastSlack=0;
int lstart = fact->nnetas + 1;
/*int nnentu = ninbas; */
int lstart_minus_nnentu=lstart-ninbas;
/* do initial column singletons - as can do faster */
for (jpivot = 1; jpivot <= nrow; ++jpivot) {
if (hincol[jpivot] == 1) {
ipivot = hrowi[mcstrt[jpivot]];
if (ipivot>lastSlack) {
lastSlack=ipivot;
} else {
/* so we can't put a structural over a slack */
break;
}
kipis = mrstrt[ipivot];
#if 1
assert (hcoli[kipis]==jpivot);
#else
if (hcoli[kipis]!=jpivot) {
kpivot=kipis+1;
while(hcoli[kpivot]!=jpivot) kpivot++;
#ifdef DEBUG
kipie = kipis + hinrow[ipivot] ;
if (kpivot>=kipie) {
abort();
}
#endif
pivot=dluval[kpivot];
dluval[kpivot] = dluval[kipis];
dluval[kipis] = pivot;
hcoli[kpivot] = hcoli[kipis];
hcoli[kipis] = jpivot;
}
#endif
if (dluval[kipis]==SLACK_VALUE) {
/* record the new pivot row and column */
++fact->npivots;
rlink[ipivot].pre = -fact->npivots;
clink[jpivot].pre = -fact->npivots;
hincol[jpivot]=0;
fact->nuspike += hinrow[ipivot];
} else {
break;
}
} else {
break;
}
}
/* Fill queue with other column singletons and clean up */
maxstk = 0;
for (j = 1; j <= nrow; ++j) {
if (hincol[j]) {
int n=0;
kcs = mcstrt[j];
kce = mcstrt[j + 1];
for (k = kcs; k < kce; ++k) {
if (! (rlink[hrowi[k]].pre < 0)) {
n++;
}
}
hincol[j] = n;
if (n == 1) {
/* we just created a new singleton column - enqueue it */
++maxstk;
stack[maxstk] = j;
}
}
}
stackc = 0; /* (1) */
while (! (stackc >= maxstk)) { /* (1) */
/* dequeue the next entry */
++stackc;
jpivot = stack[stackc];
/* (15) */
if (hincol[jpivot] != 0) {
for (k = mcstrt[jpivot]; rlink[hrowi[k]].pre < 0; k++) {
/* (4) */
}
ipivot = hrowi[k];
/* All the columns in this row are being shortened. */
kipis = mrstrt[ipivot];
kipie = kipis + hinrow[ipivot] ;
for (k = kipis; k < kipie; ++k) {
j = hcoli[k];
--hincol[j]; /* (3) (6) */
if (j == jpivot) {
kpivot = k; /* (11) */
} else if (hincol[j] == 1) {
/* we just created a new singleton column - enqueue it */
++maxstk;
stack[maxstk] = j;
}
}
/* record the new pivot row and column */
++fact->npivots;
rlink[ipivot].pre = -fact->npivots;
clink[jpivot].pre = -fact->npivots;
fact->nuspike += hinrow[ipivot];
/* check the pivot */
assert (kpivot>0);
pivot = dluval[kpivot];
if (fabs(pivot) < drtpiv) {
irtcod = 7;
++(*nsingp);
rlink[ipivot].pre = -nrow - 1;
clink[jpivot].pre = -nrow - 1;
}
/* swap the pivot column entry with the first one. */
/* I don't know why. */
dluval[kpivot] = dluval[kipis];
dluval[kipis] = pivot;
hcoli[kpivot] = hcoli[kipis];
hcoli[kipis] = jpivot;
}
}
/* (8) */
/* The entire basis may already be triangular */
if (fact->npivots < nrow) {
/* (9) */
kstart = 0;
for (j = 1; j <= nrow; ++j) {
if (! (clink[j].pre < 0)) {
kcs = mcstrt[j];
kce = mcstrt[j + 1];
mcstrt[j] = kstart + 1;
for (k = kcs; k < kce; ++k) {
if (! (rlink[hrowi[k]].pre < 0)) {
++kstart;
hrowi[kstart] = hrowi[k];
}
}
hincol[j] = kstart - mcstrt[j] + 1;
}
}
xnewco = kstart;
/* Fill stack with initial row singletons that haven't been pivoted away */
stackr = 0;
for (i = 1; i <= nrow; ++i) {
if (! (rlink[i].pre < 0) &&
(hinrow[i] == 1)) {
++stackr;
stack[stackr] = i;
}
}
while (! (stackr <= 0)) {
ipivot = stack[stackr];
assert (ipivot);
--stackr;
#if 1
assert (rlink[ipivot].pre>=0);
#else
/* This test is probably unnecessary: rlink[i].pre < 0 ==> hinrow[i]==0 */
if (rlink[ipivot].pre < 0) {
continue;
}
#endif
/* (15) */
if (hinrow[ipivot] != 0) {
/* This is a singleton row, which means it has exactly one column */
jpivot = hcoli[mrstrt[ipivot]];
kjpis = mcstrt[jpivot];
epivco = hincol[jpivot] - 1;
hincol[jpivot] = 0; /* this column is being pivoted away */
/* (11) */
kjpie = kjpis + epivco;
for (k = kjpis; k <= kjpie; ++k) {
if (ipivot == hrowi[k])
break;
}
/* ASSERT (k <= kjpie) */
/* move the last row entry for the pivot column into the pivot row's entry */
/* I don't know why */
hrowi[k] = hrowi[kjpie];
/* invalidate the (old) last row entry of the pivot column */
/* I don't know why */
hrowi[kjpie] = 0;
/* (12) */
if (! (xnewro + epivco < lstart)) {
int kstart;
if (! (epivco < lstart_minus_nnentu)) {
irtcod = -5;
break;
}
kstart = c_ekkrwco(fact,dluval, hcoli, mrstrt, hinrow, xnewro);
++ncompactions;
kmxeta += (xnewro - kstart) << 1;
xnewro = kstart;
}
if (! (xnewco + epivco < lstart)) {
if (! (epivco < lstart_minus_nnentu)) {
irtcod = -5;
break;
}
xnewco = c_ekkclco(fact,hrowi, mcstrt, hincol, xnewco);
++ncompactions;
/* HINCOL MAY HAVE CHANGED ??? (JJHF) */
epivco = hincol[jpivot];
}
/* record the new pivot row and column */
++fact->npivots;
rlink[ipivot].pre = -fact->npivots;
clink[jpivot].pre = -fact->npivots;
/* no update for nuspik */
/* check the pivot */
pivot = dluval[mrstrt[ipivot]];
if (fabs(pivot) < drtpiv) {
/* If the pivot is too small, reject it, but keep going */
irtcod = 7;
rlink[ipivot].pre = -nrow - 1;
clink[jpivot].pre = -nrow - 1;
}
/* Perform numerical part of elimination. */
if (! (epivco <= 0)) {
++xnetal;
mcstrt[xnetal] = lstart - 1;
hpivco[xnetal] = ipivot;
pivot = -1.f / pivot;
kcs = mcstrt[jpivot];
kce = kcs + epivco - 1;
hincol[jpivot] = 0;
for (kc = kcs; kc <= kce; ++kc) {
npr = hrowi[kc];
/* why bother? */
hrowi[kc] = 0;
--hinrow[npr]; /* (3) */
if (hinrow[npr] == 1) {
/* this may create new singleton rows */
++stackr;
stack[stackr] = npr;
}
/* (11) */
knprs = mrstrt[npr];
knpre = knprs + hinrow[npr];
for (k = knprs; k <= knpre; ++k) {
if (jpivot == hcoli[k]) {
kpivot = k;
break;
}
}
/* ASSERT (kpivot <= knpre) */
{
/* (10) */
double elemnt = dluval[kpivot];
dluval[kpivot] = dluval[knpre];
hcoli[kpivot] = hcoli[knpre];
hcoli[knpre] = 0; /* (14) */
/* store elementary row transformation */
--lstart;
dluval[lstart] = elemnt * pivot;
hcoli[lstart] = npr;
}
}
}
}
}
}
/* (8) */
*xnewcop = xnewco;
*xnewrop = xnewro;
fact->xnetal = xnetal;
fact->nnentu = lstart - lstart_minus_nnentu;
fact->kmxeta = kmxeta;
*ncompactionsp = ncompactions;
return (irtcod);
} /* c_ekktria */