You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
3820 lines
105 KiB
3820 lines
105 KiB
|
|
#include <string.h>
|
|
#include "commonlib.h"
|
|
#include "lp_lib.h"
|
|
#include "lp_scale.h"
|
|
#include "lp_report.h"
|
|
#include "lp_price.h"
|
|
#include "lp_pricePSE.h"
|
|
#include "lp_matrix.h"
|
|
|
|
#ifdef FORTIFY
|
|
# include "lp_fortify.h"
|
|
#endif
|
|
|
|
|
|
/* -------------------------------------------------------------------------
|
|
Basic matrix routines in lp_solve v5.0+
|
|
-------------------------------------------------------------------------
|
|
Author: Michel Berkelaar (to lp_solve v3.2),
|
|
Kjell Eikland (v4.0 and forward)
|
|
Contact: kjell.eikland@broadpark.no
|
|
License terms: LGPL.
|
|
|
|
Requires: lp_lib.h, lp_pricerPSE.h, lp_matrix.h
|
|
|
|
Release notes:
|
|
v5.0.0 1 January 2004 First integrated and repackaged version.
|
|
v5.0.1 7 May 2004 Added matrix transpose function.
|
|
v5.1.0 20 July 2004 Reworked with flexible matrix storage model.
|
|
v5.2.0 10 January 2005 Added fast deletion methods.
|
|
Added data extraction to matrix method.
|
|
Changed to explicit OF storage mode.
|
|
|
|
------------------------------------------------------------------------- */
|
|
|
|
STATIC MATrec *mat_create(lprec *lp, int rows, int columns, REAL epsvalue)
|
|
{
|
|
MATrec *newmat;
|
|
|
|
newmat = (MATrec *) calloc(1, sizeof(*newmat));
|
|
newmat->lp = lp;
|
|
|
|
newmat->rows_alloc = 0;
|
|
newmat->columns_alloc = 0;
|
|
newmat->mat_alloc = 0;
|
|
|
|
inc_matrow_space(newmat, rows);
|
|
newmat->rows = rows;
|
|
inc_matcol_space(newmat, columns);
|
|
newmat->columns = columns;
|
|
inc_mat_space(newmat, 0);
|
|
|
|
newmat->epsvalue = epsvalue;
|
|
|
|
return( newmat );
|
|
}
|
|
|
|
STATIC void mat_free(MATrec **matrix)
|
|
{
|
|
if((matrix == NULL) || (*matrix == NULL))
|
|
return;
|
|
|
|
#if MatrixColAccess==CAM_Record
|
|
FREE((*matrix)->col_mat);
|
|
#else /*if MatrixColAccess==CAM_Vector*/
|
|
FREE((*matrix)->col_mat_colnr);
|
|
FREE((*matrix)->col_mat_rownr);
|
|
FREE((*matrix)->col_mat_value);
|
|
#endif
|
|
FREE((*matrix)->col_end);
|
|
FREE((*matrix)->col_tag);
|
|
|
|
#if MatrixRowAccess==RAM_Index
|
|
FREE((*matrix)->row_mat);
|
|
#elif MatrixColAccess==CAM_Record
|
|
FREE((*matrix)->row_mat);
|
|
#else /*if MatrixRowAccess==COL_Vector*/
|
|
FREE((*matrix)->row_mat_colnr);
|
|
FREE((*matrix)->row_mat_rownr);
|
|
FREE((*matrix)->row_mat_value);
|
|
#endif
|
|
FREE((*matrix)->row_end);
|
|
FREE((*matrix)->row_tag);
|
|
|
|
FREE((*matrix)->colmax);
|
|
FREE((*matrix)->rowmax);
|
|
|
|
FREE(*matrix);
|
|
}
|
|
|
|
STATIC MYBOOL mat_memopt(MATrec *mat, int rowextra, int colextra, int nzextra)
|
|
{
|
|
MYBOOL status = TRUE;
|
|
int matalloc, colalloc, rowalloc;
|
|
|
|
if((mat == NULL) ||
|
|
#if 0
|
|
(++rowextra < 1) || (++colextra < 1) || (++nzextra < 1))
|
|
#else
|
|
(rowextra < 0) || (colextra < 0) || (nzextra < 0))
|
|
#endif
|
|
return( FALSE );
|
|
|
|
mat->rows_alloc = MIN(mat->rows_alloc, mat->rows + rowextra);
|
|
mat->columns_alloc = MIN(mat->columns_alloc, mat->columns + colextra);
|
|
mat->mat_alloc = MIN(mat->mat_alloc, mat->col_end[mat->columns] + nzextra);
|
|
#if 0
|
|
rowalloc = mat->rows_alloc;
|
|
colalloc = mat->columns_alloc;
|
|
matalloc = mat->mat_alloc;
|
|
#else
|
|
rowalloc = mat->rows_alloc + 1;
|
|
colalloc = mat->columns_alloc + 1;
|
|
matalloc = mat->mat_alloc + 1;
|
|
#endif
|
|
|
|
#if MatrixColAccess==CAM_Record
|
|
mat->col_mat = (MATitem *) realloc(mat->col_mat, matalloc * sizeof(*(mat->col_mat)));
|
|
status &= (mat->col_mat != NULL);
|
|
#else /*if MatrixColAccess==CAM_Vector*/
|
|
status &= allocINT(mat->lp, &(mat->col_mat_colnr), matalloc, AUTOMATIC) &&
|
|
allocINT(mat->lp, &(mat->col_mat_rownr), matalloc, AUTOMATIC) &&
|
|
allocREAL(mat->lp, &(mat->col_mat_value), matalloc, AUTOMATIC);
|
|
#endif
|
|
status &= allocINT(mat->lp, &mat->col_end, colalloc, AUTOMATIC);
|
|
if(mat->col_tag != NULL)
|
|
status &= allocINT(mat->lp, &mat->col_tag, colalloc, AUTOMATIC);
|
|
|
|
#if MatrixRowAccess==RAM_Index
|
|
status &= allocINT(mat->lp, &(mat->row_mat), matalloc, AUTOMATIC);
|
|
#elif MatrixColAccess==CAM_Record
|
|
mat->row_mat = (MATitem *) realloc(mat->row_mat, matalloc * sizeof(*(mat->row_mat)));
|
|
status &= (mat->row_mat != NULL);
|
|
#else /*if MatrixRowAccess==COL_Vector*/
|
|
status &= allocINT(mat->lp, &(mat->row_mat_colnr), matalloc, AUTOMATIC) &&
|
|
allocINT(mat->lp, &(mat->row_mat_rownr), matalloc, AUTOMATIC) &&
|
|
allocREAL(mat->lp, &(mat->row_mat_value), matalloc, AUTOMATIC);
|
|
#endif
|
|
status &= allocINT(mat->lp, &mat->row_end, rowalloc, AUTOMATIC);
|
|
if(mat->row_tag != NULL)
|
|
status &= allocINT(mat->lp, &mat->row_tag, rowalloc, AUTOMATIC);
|
|
|
|
if(mat->colmax != NULL)
|
|
status &= allocREAL(mat->lp, &(mat->colmax), colalloc, AUTOMATIC);
|
|
if(mat->rowmax != NULL)
|
|
status &= allocREAL(mat->lp, &(mat->rowmax), rowalloc, AUTOMATIC);
|
|
|
|
return( status );
|
|
}
|
|
|
|
STATIC MYBOOL inc_mat_space(MATrec *mat, int mindelta)
|
|
{
|
|
int spaceneeded, nz = mat_nonzeros(mat);
|
|
|
|
if(mindelta <= 0)
|
|
mindelta = MAX(mat->rows, mat->columns) + 1;
|
|
spaceneeded = DELTA_SIZE(mindelta, nz);
|
|
SETMAX(mindelta, spaceneeded);
|
|
|
|
if(mat->mat_alloc == 0)
|
|
spaceneeded = mindelta;
|
|
else
|
|
spaceneeded = nz + mindelta;
|
|
|
|
if(spaceneeded >= mat->mat_alloc) {
|
|
/* Let's allocate at least MAT_START_SIZE entries */
|
|
if(mat->mat_alloc < MAT_START_SIZE)
|
|
mat->mat_alloc = MAT_START_SIZE;
|
|
|
|
/* Increase the size by RESIZEFACTOR each time it becomes too small */
|
|
while(spaceneeded >= mat->mat_alloc)
|
|
mat->mat_alloc += mat->mat_alloc / RESIZEFACTOR;
|
|
|
|
#if MatrixColAccess==CAM_Record
|
|
mat->col_mat = (MATitem *) realloc(mat->col_mat, (mat->mat_alloc) * sizeof(*(mat->col_mat)));
|
|
#else /*if MatrixColAccess==CAM_Vector*/
|
|
allocINT(mat->lp, &(mat->col_mat_colnr), mat->mat_alloc, AUTOMATIC);
|
|
allocINT(mat->lp, &(mat->col_mat_rownr), mat->mat_alloc, AUTOMATIC);
|
|
allocREAL(mat->lp, &(mat->col_mat_value), mat->mat_alloc, AUTOMATIC);
|
|
#endif
|
|
|
|
#if MatrixRowAccess==RAM_Index
|
|
allocINT(mat->lp, &(mat->row_mat), mat->mat_alloc, AUTOMATIC);
|
|
#elif MatrixColAccess==CAM_Record
|
|
mat->row_mat = (MATitem *) realloc(mat->row_mat, (mat->mat_alloc) * sizeof(*(mat->row_mat)));
|
|
#else /*if MatrixColAccess==CAM_Vector*/
|
|
allocINT(mat->lp, &(mat->row_mat_colnr), mat->mat_alloc, AUTOMATIC);
|
|
allocINT(mat->lp, &(mat->row_mat_rownr), mat->mat_alloc, AUTOMATIC);
|
|
allocREAL(mat->lp, &(mat->row_mat_value), mat->mat_alloc, AUTOMATIC);
|
|
#endif
|
|
}
|
|
return(TRUE);
|
|
}
|
|
|
|
STATIC MYBOOL inc_matrow_space(MATrec *mat, int deltarows)
|
|
{
|
|
int rowsum, oldrowsalloc;
|
|
MYBOOL status = TRUE;
|
|
|
|
/* Adjust lp row structures */
|
|
if(mat->rows+deltarows >= mat->rows_alloc) {
|
|
|
|
/* Update memory allocation and sizes */
|
|
oldrowsalloc = mat->rows_alloc;
|
|
deltarows = DELTA_SIZE(deltarows, mat->rows);
|
|
SETMAX(deltarows, DELTAROWALLOC);
|
|
mat->rows_alloc += deltarows;
|
|
rowsum = mat->rows_alloc + 1;
|
|
|
|
/* Update row pointers */
|
|
status = allocINT(mat->lp, &mat->row_end, rowsum, AUTOMATIC);
|
|
mat->row_end_valid = FALSE;
|
|
}
|
|
return( status );
|
|
}
|
|
|
|
STATIC MYBOOL inc_matcol_space(MATrec *mat, int deltacols)
|
|
{
|
|
int i, colsum, oldcolsalloc;
|
|
MYBOOL status = TRUE;
|
|
|
|
/* Adjust lp column structures */
|
|
if(mat->columns+deltacols >= mat->columns_alloc) {
|
|
|
|
/* Update memory allocation and sizes */
|
|
oldcolsalloc = mat->columns_alloc;
|
|
deltacols = DELTA_SIZE(deltacols, mat->columns);
|
|
SETMAX(deltacols, DELTACOLALLOC);
|
|
mat->columns_alloc += deltacols;
|
|
colsum = mat->columns_alloc + 1;
|
|
status = allocINT(mat->lp, &mat->col_end, colsum, AUTOMATIC);
|
|
|
|
/* Update column pointers */
|
|
if(oldcolsalloc == 0)
|
|
mat->col_end[0] = 0;
|
|
for(i = MIN(oldcolsalloc, mat->columns) + 1; i < colsum; i++)
|
|
mat->col_end[i] = mat->col_end[i-1];
|
|
mat->row_end_valid = FALSE;
|
|
}
|
|
return( status );
|
|
}
|
|
|
|
STATIC int mat_collength(MATrec *mat, int colnr)
|
|
{
|
|
return( mat->col_end[colnr] - mat->col_end[colnr-1] );
|
|
}
|
|
|
|
STATIC int mat_rowlength(MATrec *mat, int rownr)
|
|
{
|
|
if(mat_validate(mat)) {
|
|
if(rownr <= 0)
|
|
return( mat->row_end[0] );
|
|
else
|
|
return( mat->row_end[rownr] - mat->row_end[rownr-1] );
|
|
}
|
|
else
|
|
return( 0 );
|
|
}
|
|
|
|
STATIC int mat_nonzeros(MATrec *mat)
|
|
{
|
|
return( mat->col_end[mat->columns] );
|
|
}
|
|
|
|
STATIC MYBOOL mat_indexrange(MATrec *mat, int index, MYBOOL isrow, int *startpos, int *endpos)
|
|
{
|
|
#ifdef Paranoia
|
|
if(isrow && ((index < 0) || (index > mat->rows)))
|
|
return( FALSE );
|
|
else if(!isrow && ((index < 1) || (index > mat->columns)))
|
|
return( FALSE );
|
|
#endif
|
|
|
|
if(isrow && mat_validate(mat)) {
|
|
if(index == 0)
|
|
*startpos = 0;
|
|
else
|
|
*startpos = mat->row_end[index-1];
|
|
*endpos = mat->row_end[index];
|
|
}
|
|
else {
|
|
*startpos = mat->col_end[index-1];
|
|
*endpos = mat->col_end[index];
|
|
}
|
|
return( TRUE );
|
|
}
|
|
|
|
STATIC int mat_shiftrows(MATrec *mat, int *bbase, int delta, LLrec *varmap)
|
|
{
|
|
int j, k, i, ii, thisrow, *colend, base;
|
|
MYBOOL preparecompact = FALSE;
|
|
int *rownr;
|
|
|
|
if(delta == 0)
|
|
return( 0 );
|
|
base = abs(*bbase);
|
|
|
|
if(delta > 0) {
|
|
|
|
/* Insert row by simply incrementing existing row indeces */
|
|
if(base <= mat->rows) {
|
|
k = mat_nonzeros(mat);
|
|
rownr = &COL_MAT_ROWNR(0);
|
|
for(ii = 0; ii < k; ii++, rownr += matRowColStep) {
|
|
if(*rownr >= base)
|
|
*rownr += delta;
|
|
}
|
|
}
|
|
|
|
/* Set defaults (actual basis set in separate procedure) */
|
|
for(i = 0; i < delta; i++) {
|
|
ii = base + i;
|
|
mat->row_end[ii] = 0;
|
|
}
|
|
}
|
|
else if(base <= mat->rows) {
|
|
|
|
/* Check for preparation of mass-deletion of rows */
|
|
preparecompact = (MYBOOL) (varmap != NULL);
|
|
if(preparecompact) {
|
|
/* Create the offset array */
|
|
int *newrowidx = NULL;
|
|
allocINT(mat->lp, &newrowidx, mat->rows+1, FALSE);
|
|
newrowidx[0] = 0;
|
|
delta = 0;
|
|
for(j = 1; j <= mat->rows; j++) {
|
|
if(isActiveLink(varmap, j)) {
|
|
delta++;
|
|
newrowidx[j] = delta;
|
|
}
|
|
else
|
|
newrowidx[j] = -1;
|
|
}
|
|
k = 0;
|
|
delta = 0;
|
|
base = mat_nonzeros(mat);
|
|
rownr = &COL_MAT_ROWNR(0);
|
|
for(i = 0; i < base; i++, rownr += matRowColStep) {
|
|
thisrow = newrowidx[*rownr];
|
|
if(thisrow < 0) {
|
|
*rownr = -1;
|
|
delta++;
|
|
}
|
|
else
|
|
*rownr = thisrow;
|
|
}
|
|
FREE(newrowidx);
|
|
return(delta);
|
|
}
|
|
|
|
/* Check if we should prepare for compacting later
|
|
(this is in order to speed up multiple row deletions) */
|
|
preparecompact = (MYBOOL) (*bbase < 0);
|
|
if(preparecompact)
|
|
*bbase = my_flipsign((*bbase));
|
|
|
|
/* First make sure we don't cross the row count border */
|
|
if(base-delta-1 > mat->rows)
|
|
delta = base - mat->rows - 1;
|
|
|
|
/* Then scan over all entries shifting and updating rows indeces */
|
|
if(preparecompact) {
|
|
k = 0;
|
|
for(j = 1, colend = mat->col_end + 1;
|
|
j <= mat->columns; j++, colend++) {
|
|
i = k;
|
|
k = *colend;
|
|
rownr = &COL_MAT_ROWNR(i);
|
|
for(; i < k; i++, rownr += matRowColStep) {
|
|
thisrow = *rownr;
|
|
if(thisrow < base)
|
|
continue;
|
|
else if(thisrow >= base-delta)
|
|
*rownr += delta;
|
|
else
|
|
*rownr = -1;
|
|
}
|
|
}
|
|
}
|
|
else {
|
|
k = 0;
|
|
ii = 0;
|
|
for(j = 1, colend = mat->col_end + 1;
|
|
j <= mat->columns; j++, colend++) {
|
|
i = k;
|
|
k = *colend;
|
|
rownr = &COL_MAT_ROWNR(i);
|
|
for(; i < k; i++, rownr += matRowColStep) {
|
|
thisrow = *rownr;
|
|
if(thisrow >= base) {
|
|
if(thisrow >= base-delta)
|
|
*rownr += delta;
|
|
else
|
|
continue;
|
|
}
|
|
if(ii != i) {
|
|
COL_MAT_COPY(ii, i);
|
|
}
|
|
ii++;
|
|
}
|
|
*colend = ii;
|
|
}
|
|
}
|
|
}
|
|
return( 0 );
|
|
}
|
|
|
|
/* Map-based compacting+insertion of matrix elements without changing row and column indeces.
|
|
When mat2 is NULL, a simple compacting of non-deleted rows and columns is done. */
|
|
STATIC int mat_mapreplace(MATrec *mat, LLrec *rowmap, LLrec *colmap, MATrec *mat2)
|
|
{
|
|
lprec *lp = mat->lp;
|
|
int i, ib, ie, ii, j, jj, jb, je, nz, *colend, *rownr, *rownr2, *indirect = NULL;
|
|
REAL *value, *value2;
|
|
|
|
/* Check if there is something to insert */
|
|
if((mat2 != NULL) && ((mat2->col_tag == NULL) || (mat2->col_tag[0] <= 0) || (mat_nonzeros(mat2) == 0)))
|
|
return( 0 );
|
|
|
|
/* Create map and sort by increasing index in "mat" */
|
|
if(mat2 != NULL) {
|
|
jj = mat2->col_tag[0];
|
|
allocINT(lp, &indirect, jj+1, FALSE);
|
|
indirect[0] = jj;
|
|
for(i = 1; i <= jj; i++)
|
|
indirect[i] = i;
|
|
hpsortex(mat2->col_tag, jj, 1, sizeof(*indirect), FALSE, compareINT, indirect);
|
|
}
|
|
|
|
/* Do the compacting */
|
|
mat->row_end_valid = FALSE;
|
|
nz = mat->col_end[mat->columns];
|
|
ie = 0;
|
|
ii = 0;
|
|
if((mat2 == NULL) || (indirect[0] == 0)) {
|
|
je = mat->columns + 1;
|
|
jj = 1;
|
|
jb = 0;
|
|
}
|
|
else {
|
|
je = indirect[0];
|
|
jj = 0;
|
|
do {
|
|
jj++;
|
|
jb = mat2->col_tag[jj];
|
|
} while(jb <= 0);
|
|
|
|
}
|
|
for(j = 1, colend = mat->col_end + 1;
|
|
j <= mat->columns; j++, colend++) {
|
|
ib = ie;
|
|
ie = *colend;
|
|
|
|
/* Always skip (condense) replacement columns */
|
|
if(j == jb) {
|
|
jj++;
|
|
if(jj <= je)
|
|
jb = mat2->col_tag[jj];
|
|
else
|
|
jb = mat->columns + 1;
|
|
}
|
|
|
|
/* Only include active columns */
|
|
else if(isActiveLink(colmap, j)) {
|
|
rownr = &COL_MAT_ROWNR(ib);
|
|
for(; ib < ie; ib++, rownr += matRowColStep) {
|
|
|
|
/* Also make sure the row is active */
|
|
if(isActiveLink(rowmap, *rownr)) {
|
|
if(ii != ib) {
|
|
COL_MAT_COPY(ii, ib);
|
|
}
|
|
ii++;
|
|
}
|
|
}
|
|
}
|
|
*colend = ii;
|
|
}
|
|
if(mat2 == NULL)
|
|
goto Finish;
|
|
|
|
/* Tally non-zero insertions */
|
|
i = 0;
|
|
for(j = 1; j <= mat2->col_tag[0]; j++) {
|
|
jj = mat2->col_tag[j];
|
|
if((jj > 0) && isActiveLink(colmap, jj)) {
|
|
jj = indirect[j];
|
|
je = mat2->col_end[jj];
|
|
jb = mat2->col_end[jj-1];
|
|
rownr2 = &COL_MAT2_ROWNR(jb);
|
|
for(; jb < je; jb++, rownr2 += matRowColStep) {
|
|
if((*rownr2 > 0) && isActiveLink(rowmap, *rownr2))
|
|
i++;
|
|
}
|
|
}
|
|
}
|
|
|
|
/* Make sure we have enough matrix space */
|
|
ii = mat->col_end[mat->columns] + i;
|
|
if(mat->mat_alloc <= ii)
|
|
inc_mat_space(mat, i);
|
|
|
|
/* Do shifting and insertion - loop from the end going forward */
|
|
jj = indirect[0];
|
|
jj = mat2->col_tag[jj];
|
|
for(j = mat->columns, colend = mat->col_end + mat->columns, ib = *colend;
|
|
j > 0; j--) {
|
|
|
|
/* Update indeces for this loop */
|
|
ie = ib;
|
|
*colend = ii;
|
|
colend--;
|
|
ib = *colend;
|
|
|
|
/* Insert new values */
|
|
if(j == jj) {
|
|
/* Only include an active column */
|
|
if(isActiveLink(colmap, j)) {
|
|
jj = indirect[0];
|
|
jj = indirect[jj];
|
|
rownr = &COL_MAT_ROWNR(ii-1);
|
|
value = &COL_MAT_VALUE(ii-1);
|
|
jb = mat2->col_end[jj-1];
|
|
je = mat2->col_end[jj] - 1;
|
|
rownr2 = &COL_MAT2_ROWNR(je);
|
|
value2 = &COL_MAT2_VALUE(je);
|
|
|
|
/* Process constraint coefficients */
|
|
for(; je >= jb; je--, rownr2 -= matRowColStep, value2 -= matValueStep) {
|
|
i = *rownr2;
|
|
if(i == 0) {
|
|
i = -1;
|
|
break;
|
|
}
|
|
else if(isActiveLink(rowmap, i)) {
|
|
ii--;
|
|
*rownr = i;
|
|
rownr -= matRowColStep;
|
|
*value = my_chsign(is_chsign(lp, i), *value2);
|
|
value -= matValueStep;
|
|
}
|
|
}
|
|
|
|
/* Then handle the objective */
|
|
if(i == -1) {
|
|
lp->orig_obj[j] = my_chsign(is_maxim(lp), *value2);
|
|
rownr2 -= matRowColStep;
|
|
value2 -= matValueStep;
|
|
}
|
|
else
|
|
lp->orig_obj[j] = 0;
|
|
|
|
}
|
|
/* Update replacement column index or break if no more candidates */
|
|
jj = --indirect[0];
|
|
if(jj == 0)
|
|
break;
|
|
jj = mat2->col_tag[jj];
|
|
if(jj <= 0)
|
|
break;
|
|
}
|
|
/* Shift existing values down */
|
|
else {
|
|
if(isActiveLink(colmap, j))
|
|
while(ie > ib) {
|
|
ii--;
|
|
ie--;
|
|
if(ie != ii) {
|
|
COL_MAT_COPY(ii, ie);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/* Return the delta number of non-zero elements */
|
|
Finish:
|
|
nz -= mat->col_end[mat->columns];
|
|
FREE(indirect);
|
|
|
|
return( nz );
|
|
}
|
|
|
|
/* Routines to compact rows in matrix based on precoded entries */
|
|
STATIC int mat_zerocompact(MATrec *mat)
|
|
{
|
|
return( mat_rowcompact(mat, TRUE) );
|
|
}
|
|
STATIC int mat_rowcompact(MATrec *mat, MYBOOL dozeros)
|
|
{
|
|
int i, ie, ii, j, nn, *colend, *rownr;
|
|
REAL *value;
|
|
|
|
nn = 0;
|
|
ie = 0;
|
|
ii = 0;
|
|
for(j = 1, colend = mat->col_end + 1;
|
|
j <= mat->columns; j++, colend++) {
|
|
i = ie;
|
|
ie = *colend;
|
|
rownr = &COL_MAT_ROWNR(i);
|
|
value = &COL_MAT_VALUE(i);
|
|
for(; i < ie;
|
|
i++, rownr += matRowColStep, value += matValueStep) {
|
|
if((*rownr < 0) || (dozeros && (fabs(*value) < mat->epsvalue))) {
|
|
nn++;
|
|
continue;
|
|
}
|
|
if(ii != i) {
|
|
COL_MAT_COPY(ii, i);
|
|
}
|
|
ii++;
|
|
}
|
|
*colend = ii;
|
|
}
|
|
return( nn );
|
|
}
|
|
|
|
/* Routines to compact columns and their indeces based on precoded entries */
|
|
STATIC int mat_colcompact(MATrec *mat, int prev_rows, int prev_cols)
|
|
{
|
|
int i, ii, j, k, n_del, n_sum, *colend, *newcolend, *colnr, newcolnr;
|
|
MYBOOL deleted;
|
|
lprec *lp = mat->lp;
|
|
presolveundorec *lpundo = lp->presolve_undo;
|
|
|
|
|
|
n_sum = 0;
|
|
k = 0;
|
|
ii = 0;
|
|
newcolnr = 1;
|
|
for(j = 1, colend = newcolend = mat->col_end + 1;
|
|
j <= prev_cols; j++, colend++) {
|
|
n_del = 0;
|
|
i = k;
|
|
k = *colend;
|
|
for(colnr = &COL_MAT_COLNR(i); i < k;
|
|
i++, colnr += matRowColStep) {
|
|
if(*colnr < 0) {
|
|
n_del++;
|
|
n_sum++;
|
|
continue;
|
|
}
|
|
if(ii < i) {
|
|
COL_MAT_COPY(ii, i);
|
|
}
|
|
if(newcolnr < j) {
|
|
COL_MAT_COLNR(ii) = newcolnr;
|
|
}
|
|
ii++;
|
|
}
|
|
*newcolend = ii;
|
|
|
|
deleted = (MYBOOL) (n_del > 0);
|
|
#if 1
|
|
/* Do hoops in case there was an empty column */
|
|
deleted |= (MYBOOL) (!lp->wasPresolved && (lpundo->var_to_orig[prev_rows+j] < 0));
|
|
|
|
#endif
|
|
/* Increment column variables if current column was not deleted */
|
|
if(!deleted) {
|
|
newcolend++;
|
|
newcolnr++;
|
|
}
|
|
}
|
|
return(n_sum);
|
|
}
|
|
|
|
STATIC int mat_shiftcols(MATrec *mat, int *bbase, int delta, LLrec *varmap)
|
|
{
|
|
int i, ii, k, n, base;
|
|
|
|
|
|
k = 0;
|
|
if(delta == 0)
|
|
return( k );
|
|
base = abs(*bbase);
|
|
|
|
if(delta > 0) {
|
|
/* Shift pointers right */
|
|
for(ii = mat->columns; ii > base; ii--) {
|
|
i = ii + delta;
|
|
mat->col_end[i] = mat->col_end[ii];
|
|
}
|
|
/* Set defaults */
|
|
for(i = 0; i < delta; i++) {
|
|
ii = base + i;
|
|
mat->col_end[ii] = mat->col_end[ii-1];
|
|
}
|
|
}
|
|
else {
|
|
|
|
/* Check for preparation of mass-deletion of columns */
|
|
MYBOOL preparecompact = (MYBOOL) (varmap != NULL);
|
|
if(preparecompact) {
|
|
/* Create the offset array */
|
|
int j, *colnr, *colend;
|
|
n = 0;
|
|
k = 0;
|
|
base = 0;
|
|
for(j = 1, colend = mat->col_end + 1;
|
|
j <= mat->columns; j++, colend++) {
|
|
i = k;
|
|
k = *colend;
|
|
if(isActiveLink(varmap, j)) {
|
|
base++;
|
|
ii = base;
|
|
}
|
|
else
|
|
ii = -1;
|
|
if(ii < 0)
|
|
n += k - i;
|
|
colnr = &COL_MAT_COLNR(i);
|
|
for(; i < k; i++, colnr += matRowColStep)
|
|
*colnr = ii;
|
|
}
|
|
return(n);
|
|
}
|
|
|
|
/* Check if we should prepare for compacting later
|
|
(this is in order to speed up multiple column deletions) */
|
|
preparecompact = (MYBOOL) (*bbase < 0);
|
|
if(preparecompact)
|
|
*bbase = my_flipsign((*bbase));
|
|
|
|
/* First make sure we don't cross the column count border */
|
|
if(base-delta-1 > mat->columns)
|
|
delta = base - mat->columns - 1;
|
|
|
|
/* Then scan over all entries shifting and updating column indeces */
|
|
if(preparecompact) {
|
|
int *colnr;
|
|
n = 0;
|
|
i = mat->col_end[base-1];
|
|
k = mat->col_end[base-delta-1];
|
|
for(colnr = &COL_MAT_COLNR(i); i < k;
|
|
i++, colnr += matRowColStep) {
|
|
n++;
|
|
*colnr = -1;
|
|
}
|
|
k = n;
|
|
}
|
|
else {
|
|
/* Delete sparse matrix data, if required */
|
|
if(base <= mat->columns) {
|
|
|
|
i = mat->col_end[base-1]; /* Beginning of data to be deleted */
|
|
ii = mat->col_end[base-delta-1]; /* Beginning of data to be shifted left */
|
|
n = mat_nonzeros(mat); /* Total number of non-zeros */
|
|
k = ii-i; /* Number of entries to be deleted */
|
|
if((k > 0) && (n > i)) {
|
|
n -= ii;
|
|
COL_MAT_MOVE(i, ii, n);
|
|
}
|
|
|
|
/* Update indexes */
|
|
for(i = base; i <= mat->columns + delta; i++) {
|
|
ii = i - delta;
|
|
mat->col_end[i] = mat->col_end[ii] - k;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return( k );
|
|
}
|
|
|
|
STATIC MATrec *mat_extractmat(MATrec *mat, LLrec *rowmap, LLrec *colmap, MYBOOL negated)
|
|
{
|
|
int *rownr, *colnr, xa, na;
|
|
REAL *value;
|
|
MATrec *newmat = mat_create(mat->lp, mat->rows, mat->columns, mat->epsvalue);
|
|
|
|
/* Initialize */
|
|
na = mat_nonzeros(mat);
|
|
rownr = &COL_MAT_ROWNR(0);
|
|
colnr = &COL_MAT_COLNR(0);
|
|
value = &COL_MAT_VALUE(0);
|
|
|
|
/* Loop over the indeces, picking out values in qualifying rows and colums
|
|
(note that the loop could be speeded up for dense matrices by making an
|
|
outer loop for columns and inner loop for rows) */
|
|
for(xa = 0; xa < na; xa++, rownr += matRowColStep, colnr += matRowColStep, value += matValueStep) {
|
|
if((isActiveLink(colmap, *colnr) ^ negated) &&
|
|
(isActiveLink(rowmap, *rownr) ^ negated))
|
|
mat_setvalue(newmat, *rownr, *colnr, *value, FALSE);
|
|
}
|
|
|
|
/* Return the populated new matrix */
|
|
return( newmat );
|
|
}
|
|
|
|
STATIC MYBOOL mat_setcol(MATrec *mat, int colno, int count, REAL *column, int *rowno, MYBOOL doscale, MYBOOL checkrowmode)
|
|
{
|
|
int i, jj = 0, elmnr, orignr, newnr, firstrow;
|
|
MYBOOL *addto = NULL, isA, isNZ;
|
|
REAL value, saved = 0;
|
|
lprec *lp = mat->lp;
|
|
|
|
/* Check if we are in row order mode and should add as row instead;
|
|
the matrix will be transposed at a later stage */
|
|
if(checkrowmode && mat->is_roworder)
|
|
return( mat_setrow(mat, colno, count, column, rowno, doscale, FALSE) );
|
|
|
|
/* Initialize and validate */
|
|
isA = (MYBOOL) (mat == mat->lp->matA);
|
|
isNZ = (MYBOOL) (rowno != NULL);
|
|
if(!isNZ)
|
|
count = mat->lp->rows;
|
|
else if((count < 0) || (count > mat->rows+((mat->is_roworder) ? 0 : 1)))
|
|
return( FALSE );
|
|
if(isNZ && (count > 0)) {
|
|
if(count > 1)
|
|
sortREALByINT(column, rowno, count, 0, TRUE);
|
|
if((rowno[0] < 0) || (rowno[count-1] > mat->rows))
|
|
return( FALSE );
|
|
}
|
|
|
|
/* Capture OF definition in column mode */
|
|
if(isA && !mat->is_roworder) {
|
|
if(isNZ && (count > 0) && (rowno[0] == 0)) {
|
|
value = column[0];
|
|
#ifdef DoMatrixRounding
|
|
value = roundToPrecision(value, mat->epsvalue);
|
|
#endif
|
|
if(doscale)
|
|
value = scaled_mat(lp, value, 0, colno);
|
|
value = my_chsign(is_maxim(lp), value);
|
|
lp->orig_obj[colno] = value;
|
|
count--;
|
|
column++;
|
|
rowno++;
|
|
}
|
|
else if(!isNZ && (column[0] != 0)) {
|
|
value = saved = column[0];
|
|
#ifdef DoMatrixRounding
|
|
value = roundToPrecision(value, mat->epsvalue);
|
|
#endif
|
|
if(doscale)
|
|
value = scaled_mat(lp, value, 0, colno);
|
|
value = my_chsign(is_maxim(lp), value);
|
|
lp->orig_obj[colno] = value;
|
|
column[0] = 0;
|
|
}
|
|
else
|
|
lp->orig_obj[colno] = 0;
|
|
}
|
|
|
|
/* Optionally tally and map the new non-zero values */
|
|
firstrow = mat->rows + 1;
|
|
if(isNZ) {
|
|
newnr = count;
|
|
if(newnr) {
|
|
firstrow = rowno[0];
|
|
jj = rowno[newnr - 1];
|
|
}
|
|
}
|
|
else {
|
|
newnr = 0;
|
|
if(!allocMYBOOL(lp, &addto, mat->rows + 1, TRUE)) {
|
|
return( FALSE );
|
|
}
|
|
for(i = mat->rows; i >= 0; i--) {
|
|
if(fabs(column[i]) > mat->epsvalue) {
|
|
addto[i] = TRUE;
|
|
firstrow = i;
|
|
newnr++;
|
|
}
|
|
}
|
|
}
|
|
|
|
/* Make sure we have enough matrix space */
|
|
if(!inc_mat_space(mat, newnr)) {
|
|
newnr = 0;
|
|
goto Done;
|
|
}
|
|
|
|
/* Shift existing column data and adjust position indeces */
|
|
orignr = mat_collength(mat, colno);
|
|
elmnr = newnr - orignr;
|
|
i = mat_nonzeros(mat) - mat->col_end[colno];
|
|
if((elmnr != 0) && (i > 0)) {
|
|
COL_MAT_MOVE(mat->col_end[colno] + elmnr, mat->col_end[colno], i);
|
|
}
|
|
if(elmnr != 0)
|
|
for(i = colno; i <= mat->columns; i++)
|
|
mat->col_end[i] += elmnr;
|
|
|
|
/* We are now ready to copy the new data */
|
|
jj = mat->col_end[colno-1];
|
|
if(isNZ) {
|
|
for(i = 0; i < count; jj++, i++) {
|
|
value = column[i];
|
|
#ifdef DoMatrixRounding
|
|
value = roundToPrecision(value, mat->epsvalue);
|
|
#endif
|
|
if(mat->is_roworder) { /* Fix following Ingmar Stein bug report 12.10.2006 */
|
|
if(isA && doscale)
|
|
value = scaled_mat(lp, value, colno, rowno[i]);
|
|
if(isA)
|
|
value = my_chsign(is_chsign(lp, colno), value);
|
|
}
|
|
else {
|
|
if(isA && doscale)
|
|
value = scaled_mat(lp, value, rowno[i], colno);
|
|
if(isA)
|
|
value = my_chsign(is_chsign(lp, rowno[i]), value);
|
|
}
|
|
SET_MAT_ijA(jj, rowno[i], colno, value);
|
|
}
|
|
}
|
|
else {
|
|
for(i = firstrow; i <= mat->rows; i++) {
|
|
if(!addto[i])
|
|
continue;
|
|
value = column[i];
|
|
#ifdef DoMatrixRounding
|
|
value = roundToPrecision(value, mat->epsvalue);
|
|
#endif
|
|
if(mat->is_roworder) { /* Fix following Ingmar Stein bug report 12.10.2006 */
|
|
if(isA && doscale)
|
|
value = scaled_mat(lp, value, colno, i);
|
|
if(isA)
|
|
value = my_chsign(is_chsign(lp, colno), value);
|
|
}
|
|
else {
|
|
if(isA && doscale)
|
|
value = scaled_mat(lp, value, i, colno);
|
|
if(isA)
|
|
value = my_chsign(is_chsign(lp, i), value);
|
|
}
|
|
SET_MAT_ijA(jj, i, colno, value);
|
|
jj++;
|
|
}
|
|
}
|
|
mat->row_end_valid = FALSE;
|
|
|
|
/* Finish and return */
|
|
Done:
|
|
if(saved != 0)
|
|
column[0] = saved;
|
|
FREE(addto);
|
|
return( TRUE );
|
|
|
|
}
|
|
|
|
STATIC MYBOOL mat_mergemat(MATrec *target, MATrec *source, MYBOOL usecolmap)
|
|
{
|
|
lprec *lp = target->lp;
|
|
int i, ix, iy, n, *colmap = NULL;
|
|
REAL *colvalue = NULL;
|
|
|
|
if((target->rows < source->rows) || !allocREAL(lp, &colvalue, target->rows+1, FALSE))
|
|
return( FALSE );
|
|
|
|
if(usecolmap) {
|
|
n = source->col_tag[0];
|
|
allocINT(lp, &colmap, n+1, FALSE);
|
|
for(i = 1; i <= n; i++)
|
|
colmap[i] = i;
|
|
hpsortex(source->col_tag, n, 1, sizeof(*colmap), FALSE, compareINT, colmap);
|
|
}
|
|
else
|
|
n = source->columns;
|
|
for(i = 1; i <= n; i++) {
|
|
if(!usecolmap && (mat_collength(source, i) == 0))
|
|
continue;
|
|
if(usecolmap) {
|
|
ix = colmap[i];
|
|
if(ix <= 0)
|
|
continue;
|
|
iy = source->col_tag[i];
|
|
if(iy <= 0)
|
|
continue;
|
|
}
|
|
else
|
|
ix = iy = i;
|
|
mat_expandcolumn(source, ix, colvalue, NULL, FALSE);
|
|
mat_setcol(target, iy, 0, colvalue, NULL, FALSE, FALSE);
|
|
}
|
|
|
|
FREE( colvalue );
|
|
FREE( colmap );
|
|
|
|
return( TRUE );
|
|
}
|
|
|
|
STATIC int mat_nz_unused(MATrec *mat)
|
|
{
|
|
return( mat->mat_alloc - mat->col_end[mat->columns] );
|
|
}
|
|
|
|
#if 0
|
|
STATIC MYBOOL mat_setrow(MATrec *mat, int rowno, int count, REAL *row, int *colno, MYBOOL doscale, MYBOOL checkrowmode)
|
|
{
|
|
lprec *lp = mat->lp;
|
|
int delta;
|
|
int k, kk, i, ii, j, jj = 0, jj_j, elmnr, orignr, newnr, firstcol, rownr, colnr, matz = 0;
|
|
MYBOOL *addto = NULL, isA, isNZ;
|
|
REAL value = 0.0, saved = 0;
|
|
|
|
/* Check if we are in row order mode and should add as column instead;
|
|
the matrix will be transposed at a later stage */
|
|
if(checkrowmode && mat->is_roworder)
|
|
return( mat_setcol(mat, rowno, count, row, colno, doscale, FALSE) );
|
|
|
|
/* Do initialization and validation */
|
|
if(!mat_validate(mat))
|
|
return( FALSE );
|
|
isA = (MYBOOL) (mat == lp->matA);
|
|
isNZ = (MYBOOL) (colno != NULL);
|
|
if(!isNZ)
|
|
count = mat->columns;
|
|
else if((count < 0) || (count > mat->columns))
|
|
return( FALSE );
|
|
if(isNZ && (count > 0)) {
|
|
if(count > 1)
|
|
sortREALByINT(row, (int *) colno, count, 0, TRUE);
|
|
if((colno[0] < 1) || (colno[count-1] > mat->columns))
|
|
return( FALSE );
|
|
}
|
|
|
|
/* Capture OF definition in row mode */
|
|
if(isA && mat->is_roworder) {
|
|
lp->orig_obj[rowno] = 0;
|
|
if(isNZ && (count > 0) && (colno[0] == 0)) {
|
|
value = row[0];
|
|
if(doscale)
|
|
value = scaled_mat(lp, value, 0, rowno);
|
|
value = my_chsign(is_maxim(lp), value);
|
|
#ifdef DoMatrixRounding
|
|
value = roundToPrecision(value, mat->epsvalue);
|
|
#endif
|
|
lp->orig_obj[rowno] = value;
|
|
count--;
|
|
row++;
|
|
colno++;
|
|
}
|
|
else if(!isNZ && (row[0] != 0)) {
|
|
value = saved = row[0];
|
|
if(doscale)
|
|
value = scaled_mat(lp, value, 0, rowno);
|
|
value = my_chsign(is_maxim(lp), value);
|
|
#ifdef DoMatrixRounding
|
|
value = roundToPrecision(value, mat->epsvalue);
|
|
#endif
|
|
lp->orig_obj[rowno] = value;
|
|
row[0] = 0;
|
|
}
|
|
else {
|
|
lp->orig_obj[rowno] = 0;
|
|
value = 0;
|
|
}
|
|
}
|
|
|
|
/* Optionally tally and map the new non-zero values */
|
|
i = mat->row_end[rowno-1];
|
|
ii = mat->row_end[rowno]; // ****** KE 20070106 - was "-1"
|
|
firstcol = mat->columns + 1;
|
|
if(isNZ) {
|
|
/* See if we can do fast in-place replacements of leading items */
|
|
colnr = 1; /* initialise in case of an empty row */
|
|
while((i < ii) /* && (count > 0) */ && ((colnr = ROW_MAT_COLNR(i)) == *colno) && (count > 0)) {
|
|
value = *row; // ****** KE 20080111 - Added line
|
|
if(mat->is_roworder) {
|
|
if(isA && doscale)
|
|
value = scaled_mat(lp, value, colnr, rowno);
|
|
if(isA)
|
|
value = my_chsign(is_chsign(lp, colnr), value);
|
|
}
|
|
else {
|
|
if(isA && doscale)
|
|
value = scaled_mat(lp, value, rowno, colnr);
|
|
if(isA)
|
|
value = my_chsign(is_chsign(lp, rowno), value);
|
|
}
|
|
#ifdef DoMatrixRounding
|
|
value = roundToPrecision(value, mat->epsvalue);
|
|
if(value == 0)
|
|
matz++;
|
|
#endif
|
|
ROW_MAT_VALUE(i) = value;
|
|
i++;
|
|
count--;
|
|
row++;
|
|
colno++;
|
|
}
|
|
if(i >= ii)
|
|
colnr = 0;
|
|
/* Proceed with remaining entries */
|
|
newnr = count;
|
|
if(newnr > 0)
|
|
firstcol = colno[0];
|
|
}
|
|
else {
|
|
newnr = 0;
|
|
kk = mat->columns;
|
|
if(i < ii)
|
|
colnr = ROW_MAT_COLNR(i);
|
|
else
|
|
colnr = 0;
|
|
for(k = 1; k <= kk; k++) {
|
|
value = row[k]; // ****** KE 20080111 - Added line
|
|
if(fabs(value) > mat->epsvalue) {
|
|
/* See if we can do fast in-place replacements of leading items */
|
|
if((addto == NULL) && (i < ii) && (colnr == k)) {
|
|
if(mat->is_roworder) {
|
|
if(isA && doscale)
|
|
value = scaled_mat(lp, value, colnr, rowno);
|
|
if(isA)
|
|
value = my_chsign(is_chsign(lp, colnr), value);
|
|
}
|
|
else {
|
|
if(isA && doscale)
|
|
value = scaled_mat(lp, value, rowno, colnr);
|
|
if(isA)
|
|
value = my_chsign(is_chsign(lp, rowno), value);
|
|
}
|
|
#ifdef DoMatrixRounding
|
|
value = roundToPrecision(value, mat->epsvalue);
|
|
if(value == 0)
|
|
matz++;
|
|
#endif
|
|
ROW_MAT_VALUE(i) = value;
|
|
i++;
|
|
if(i < ii)
|
|
colnr = ROW_MAT_COLNR(i);
|
|
else
|
|
colnr = 0;
|
|
}
|
|
/* Otherwise update addto-list */
|
|
else {
|
|
if(addto == NULL) {
|
|
if(!allocMYBOOL(lp, &addto, mat->columns + 1, TRUE))
|
|
return( FALSE );
|
|
firstcol = k;
|
|
}
|
|
addto[k] = TRUE;
|
|
newnr++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
if(newnr == 0)
|
|
if (FALSE)
|
|
return( TRUE );
|
|
|
|
/* Make sure we have enough matrix space */
|
|
if((newnr > 0) && (mat_nz_unused(mat) <= newnr) && !inc_mat_space(mat, newnr)) {
|
|
newnr = 0;
|
|
goto Done;
|
|
}
|
|
|
|
/* Pack initial entries if existing row data has a lower column
|
|
start index than the first index of the new vector */
|
|
orignr = mat_nonzeros(mat);
|
|
/* delta = newnr - mat_rowlength(mat, rowno);*/
|
|
kk = 0;
|
|
if(rowno == 0)
|
|
ii = 0;
|
|
else
|
|
ii = mat->row_end[rowno-1];
|
|
|
|
if((orignr == 0) || (ii >= orignr))
|
|
j = firstcol;
|
|
else if(isNZ||TRUE)
|
|
j = colnr;
|
|
else
|
|
j = ROW_MAT_COLNR(ii); /* first column with a value on that row */
|
|
|
|
jj = mat->col_end[firstcol-1]; /* Set the index of the insertion point for the first new value */
|
|
if(jj >= orignr)
|
|
colnr = firstcol;
|
|
else
|
|
colnr = COL_MAT_COLNR(jj); /* first column with a value starting from firstcol */
|
|
|
|
if((j > 0) && (j < colnr)) {
|
|
jj = elmnr = mat->col_end[j-1];
|
|
for( ; j < colnr; j++) {
|
|
/* Shift entries in current column */
|
|
k = mat->col_end[j];
|
|
for( ; jj < k; jj++) {
|
|
if(COL_MAT_ROWNR(jj) != rowno) {
|
|
COL_MAT_COPY(elmnr, jj);
|
|
elmnr++;
|
|
}
|
|
}
|
|
/* Update next column start index */
|
|
mat->col_end[j] = elmnr;
|
|
}
|
|
delta = elmnr - jj; /* The shrinkage count */
|
|
}
|
|
else {
|
|
delta = 0;
|
|
/* Adjust for case where we simply append values - jj is initially the first column item */
|
|
if((mat->col_end[firstcol] == orignr) && 0)
|
|
jj = orignr;
|
|
}
|
|
|
|
/* Make sure we have sufficient space for any additional entries and move existing data down;
|
|
this ensures that we only have to relocate matrix elements up in the next stage */
|
|
jj_j = MAX(0, newnr + delta);
|
|
if(jj_j > 0) {
|
|
if(!inc_mat_space(mat, jj_j)) {
|
|
FREE(addto);
|
|
return( FALSE );
|
|
}
|
|
if(orignr-jj > 0) {
|
|
COL_MAT_MOVE(jj+jj_j, jj, orignr-jj);
|
|
}
|
|
jj += jj_j;
|
|
}
|
|
|
|
/* Handle case where the matrix was empty before (or we can simply append) */
|
|
if((delta >= 0) && (mat->col_end[firstcol] == orignr) && 0) {
|
|
if(isNZ)
|
|
elmnr = count;
|
|
else
|
|
elmnr = mat->columns;
|
|
jj_j = mat->col_end[firstcol];
|
|
for(newnr = 0; newnr < elmnr; newnr++) {
|
|
if(isNZ)
|
|
colnr = colno[newnr];
|
|
else
|
|
colnr = newnr + 1;
|
|
/* Update column start position if we have crossed a column */
|
|
while(colnr > firstcol) {
|
|
mat->col_end[firstcol] = jj_j;
|
|
firstcol++;
|
|
}
|
|
if(isNZ || ((addto != NULL) && addto[colnr])) {
|
|
if(isNZ)
|
|
value = row[newnr];
|
|
else
|
|
value = row[colnr];
|
|
if(isA && doscale)
|
|
value = scaled_mat(lp, value, rowno, colnr);
|
|
if(isA)
|
|
value = my_chsign(is_chsign(lp, rowno), value);
|
|
#ifdef DoMatrixRounding
|
|
value = roundToPrecision(value, mat->epsvalue);
|
|
if(value == 0)
|
|
matz++;
|
|
#endif
|
|
SET_MAT_ijA(jj_j, rowno, colnr, value);
|
|
jj_j++;
|
|
/* Update last column start position */
|
|
mat->col_end[firstcol] = jj_j;
|
|
firstcol++;
|
|
}
|
|
}
|
|
|
|
/* Make sure we update tail empty column offsets */
|
|
while(firstcol <= mat->columns) {
|
|
mat->col_end[firstcol] = jj_j;
|
|
firstcol++;
|
|
}
|
|
jj_j = 0;
|
|
}
|
|
|
|
/* Start from the top of the first non-zero column of the new row */
|
|
elmnr = orignr + jj_j;
|
|
if(jj <= elmnr) {
|
|
if(isNZ)
|
|
newnr = 0;
|
|
else
|
|
newnr = firstcol - 1;
|
|
j = jj - mat->col_end[firstcol-1];
|
|
colnr = firstcol;
|
|
while((jj < elmnr) || (newnr < count)) {
|
|
|
|
/* Update column start position if we have crossed a column */
|
|
while(colnr > firstcol) {
|
|
mat->col_end[firstcol] = kk;
|
|
firstcol++;
|
|
}
|
|
|
|
/* See if we have a row equal to or greater than the target row */
|
|
jj_j = jj - j;
|
|
if(jj < elmnr) {
|
|
rownr = COL_MAT_ROWNR(jj);
|
|
colnr = COL_MAT_COLNR(jj);
|
|
}
|
|
else {
|
|
rownr = rowno;
|
|
if(!isNZ) /* KE added this conditional on 13.9.2006 */
|
|
colnr = firstcol + 1;
|
|
else
|
|
colnr = mat->columns + 1;
|
|
}
|
|
|
|
if(isNZ) {
|
|
if(newnr < count)
|
|
kk = colno[newnr];
|
|
else
|
|
kk = mat->columns + 1;
|
|
}
|
|
else
|
|
kk = newnr + 1;
|
|
|
|
/* Test if there is an available new item ... */
|
|
if((isNZ && (kk > colnr)) || /* If this is not the case */
|
|
(!isNZ && ((kk > colnr) || (!addto[kk])))) {
|
|
/* DELETE if there is an existing value */
|
|
if(!isNZ && (kk <= colnr))
|
|
newnr++;
|
|
if(rownr == rowno) {
|
|
kk = jj_j;
|
|
j++;
|
|
jj++;
|
|
continue;
|
|
}
|
|
/* KEEP otherwise and move entry up */
|
|
if(!isNZ && (colnr > kk)) {
|
|
colnr = kk;
|
|
kk = jj_j;
|
|
continue;
|
|
}
|
|
}
|
|
else if((colnr > kk) || /* Existing column index > new => INSERT */
|
|
((colnr == kk) && (rownr >= rowno)) ) { /* Same column index, existing row >= target row => INSERT/REPLACE */
|
|
|
|
if(isNZ)
|
|
value = row[newnr];
|
|
else
|
|
value = row[newnr+1];
|
|
newnr++;
|
|
if(isA && doscale)
|
|
value = scaled_mat(lp, value, rowno, kk);
|
|
if(isA)
|
|
value = my_chsign(is_chsign(lp, rowno), value);
|
|
#ifdef DoMatrixRounding
|
|
value = roundToPrecision(value, mat->epsvalue);
|
|
if(value == 0)
|
|
matz++;
|
|
#endif
|
|
SET_MAT_ijA(jj_j, rowno, kk, value);
|
|
|
|
/* Adjust if we have inserted an element */
|
|
if((colnr > kk) || (rownr > rowno)) {
|
|
j--;
|
|
jj--;
|
|
}
|
|
colnr = kk;
|
|
kk = jj_j;
|
|
jj++;
|
|
continue;
|
|
}
|
|
|
|
/* Shift the matrix element up by the active difference */
|
|
if(jj_j != jj) {
|
|
COL_MAT_COPY(jj_j, jj);
|
|
}
|
|
kk = jj_j;
|
|
jj++;
|
|
|
|
}
|
|
|
|
/* Update pending / incomplete column start position */
|
|
while(colnr > firstcol) {
|
|
mat->col_end[firstcol] = kk;
|
|
firstcol++;
|
|
}
|
|
|
|
/* Make sure we update tail column offsets */
|
|
jj_j = jj - j;
|
|
while(firstcol <= mat->columns) {
|
|
mat->col_end[firstcol] = jj_j;
|
|
firstcol++;
|
|
}
|
|
}
|
|
|
|
/* Compact in the case that we added zeros and set flag for row index update */
|
|
if(matz > 0)
|
|
mat_zerocompact(mat);
|
|
mat->row_end_valid = FALSE;
|
|
|
|
Done:
|
|
if(saved != 0)
|
|
row[0] = saved;
|
|
FREE(addto);
|
|
return( (MYBOOL) (newnr > 0) );
|
|
|
|
}
|
|
|
|
#else
|
|
|
|
STATIC MYBOOL mat_setrow(MATrec *mat, int rowno, int count, REAL *row, int *colno, MYBOOL doscale, MYBOOL checkrowmode)
|
|
{
|
|
lprec *lp = mat->lp;
|
|
int delta, delta1;
|
|
int k, i, ii, j, jj_j, lendense,
|
|
origidx = 0, newidx, orignz, newnz,
|
|
rownr, colnr, colnr1;
|
|
MYBOOL isA, isNZ;
|
|
REAL value = 0.0;
|
|
|
|
/* Check if we are in row order mode and should add as column instead;
|
|
the matrix will be transposed at a later stage */
|
|
if(checkrowmode && mat->is_roworder)
|
|
return( mat_setcol(mat, rowno, count, row, colno, doscale, FALSE) );
|
|
|
|
/* Do initialization and validation */
|
|
if(!mat_validate(mat))
|
|
return( FALSE );
|
|
isA = (MYBOOL) (mat == lp->matA);
|
|
if(doscale && isA && !lp->scaling_used)
|
|
doscale = FALSE;
|
|
isNZ = (MYBOOL) (colno != NULL);
|
|
lendense = (mat->is_roworder ? lp->rows : lp->columns);
|
|
if((count < 0) || (count > lendense))
|
|
return( FALSE );
|
|
colnr1 = lendense + 1;
|
|
|
|
/* Capture OF definition in row mode */
|
|
if(isA && mat->is_roworder) {
|
|
lp->orig_obj[rowno] = 0;
|
|
if((count > 0) && (colno[0] == 0)) {
|
|
value = row[0];
|
|
if(doscale)
|
|
value = scaled_mat(lp, value, 0, rowno);
|
|
value = my_chsign(is_maxim(lp), value);
|
|
#ifdef DoMatrixRounding
|
|
value = roundToPrecision(value, mat->epsvalue);
|
|
#endif
|
|
lp->orig_obj[rowno] = value;
|
|
if(isNZ) {
|
|
colno++;
|
|
row++;
|
|
count--;
|
|
}
|
|
}
|
|
else {
|
|
lp->orig_obj[rowno] = 0;
|
|
value = 0;
|
|
}
|
|
}
|
|
|
|
/* Make local working data copies */
|
|
if(!isNZ) {
|
|
REAL *tmprow = NULL;
|
|
if(!allocINT(lp, &colno, lendense+1, FALSE))
|
|
return( FALSE );
|
|
newnz = 0;
|
|
for(i = 1; i <= lendense; i++)
|
|
if((value = row[i]) != 0) {
|
|
if((tmprow == NULL) && !allocREAL(lp, &tmprow, lendense-i+1, FALSE)) {
|
|
FREE(colno);
|
|
return( FALSE );
|
|
}
|
|
tmprow[newnz] = value;
|
|
colno[newnz++] = i;
|
|
}
|
|
count = newnz;
|
|
row = tmprow;
|
|
}
|
|
else {
|
|
int *tmpcolno = NULL;
|
|
if(!allocINT(lp, &tmpcolno, lendense, FALSE))
|
|
return( FALSE );
|
|
newnz = count;
|
|
MEMCOPY(tmpcolno, colno, newnz);
|
|
colno = tmpcolno;
|
|
if(newnz > 1)
|
|
sortREALByINT(row, (int *) colno, newnz, 0, TRUE);
|
|
if((newnz > 0) && ((colno[0] < 0) || (colno[newnz-1] > lendense))) {
|
|
FREE(colno);
|
|
newnz = 0;
|
|
return( FALSE );
|
|
}
|
|
}
|
|
|
|
/* Make sure we have enough matrix space */
|
|
i = mat->row_end[rowno-1];
|
|
ii = mat->row_end[rowno];
|
|
delta1 = delta = count - (ii-i);
|
|
colnr1 = (newnz > 0 ? colno[0] : lendense+1);
|
|
|
|
/* Pack initial entries if existing row data has a lower column
|
|
start index than the first index of the new vector */
|
|
orignz = mat_nonzeros(mat);
|
|
j = (i >= orignz ? colnr1 : ROW_MAT_COLNR(i));
|
|
|
|
/* Index of the column-top insertion point for the first new value */
|
|
origidx = mat->col_end[colnr1-1];
|
|
colnr = (origidx >= orignz ? colnr1 : COL_MAT_COLNR(origidx));
|
|
|
|
if(j < colnr) {
|
|
origidx = newidx = mat->col_end[j-1];
|
|
for( ; j < colnr; j++) {
|
|
/* Shift entries in current column */
|
|
jj_j = mat->col_end[j];
|
|
for( ; origidx < jj_j; origidx++) {
|
|
if(COL_MAT_ROWNR(origidx) != rowno) {
|
|
if(newidx != origidx) {
|
|
COL_MAT_COPY(newidx, origidx);
|
|
}
|
|
newidx++;
|
|
}
|
|
}
|
|
/* Update next column start index */
|
|
mat->col_end[j] = newidx;
|
|
}
|
|
delta = newidx - origidx; /* The first stage element shrinkage count */
|
|
}
|
|
else {
|
|
delta = 0;
|
|
newidx = origidx;
|
|
}
|
|
|
|
/* Make sure we have sufficient space for any additional entries and move existing data down;
|
|
this ensures that we only have to relocate matrix elements up in the next stage */
|
|
jj_j = MAX(0, (int) newnz + delta);
|
|
|
|
j = !((orignz == lendense) && (newnz == orignz) && (delta1 == 0)) && (jj_j > 0) && (orignz > origidx);
|
|
|
|
if ((j) && (jj_j > delta1))
|
|
delta1 = jj_j;
|
|
|
|
if((delta1 > 0) && (mat_nz_unused(mat) <= delta1) && !inc_mat_space(mat, delta1)) {
|
|
newnz = 0;
|
|
goto Done;
|
|
}
|
|
|
|
if(j) {
|
|
COL_MAT_MOVE(origidx+jj_j, origidx, orignz-origidx);
|
|
origidx += jj_j;
|
|
orignz += jj_j;
|
|
}
|
|
|
|
/* Start from the top of the first non-zero column of the new row */
|
|
newnz = 0;
|
|
j = origidx - mat->col_end[colnr1-1];
|
|
k = colnr1; /* Last column for which col_end is valid/updated */
|
|
while((colnr1 <= lendense) || (origidx < orignz)) {
|
|
|
|
/* Get the column index of the active update item */
|
|
if(newnz < count)
|
|
colnr1 = colno[newnz];
|
|
else
|
|
colnr1 = lendense + 1;
|
|
|
|
/* Get coordinate of active existing matrix entries */
|
|
if(origidx < orignz) {
|
|
rownr = COL_MAT_ROWNR(origidx);
|
|
colnr = COL_MAT_COLNR(origidx);
|
|
}
|
|
else {
|
|
if(colnr1 > lendense)
|
|
break;
|
|
rownr = rowno;
|
|
colnr = lendense + 1;
|
|
}
|
|
|
|
/* Update column start position if we just crossed into a column */
|
|
jj_j = origidx - j;
|
|
i = MIN(colnr, colnr1);
|
|
for(; k < i; k++)
|
|
mat->col_end[k] = jj_j;
|
|
|
|
/* Test if there is an available new item ... */
|
|
if(colnr1 > colnr) { /* If this is not the case */
|
|
/* DELETE if there is an existing value */
|
|
if(rownr == rowno) {
|
|
ForceDelete:
|
|
j++;
|
|
delta--;
|
|
origidx++;
|
|
continue;
|
|
}
|
|
}
|
|
else if((colnr > colnr1) || /* Existing column index > new => INSERT */
|
|
((colnr == colnr1) && (rownr >= rowno)) ) { /* Same column index, existing row >= target row => INSERT/REPLACE */
|
|
|
|
value = row[newnz];
|
|
newnz++;
|
|
if(isA && doscale)
|
|
value = scaled_mat(lp, value, rowno, colnr1);
|
|
if(isA)
|
|
value = my_chsign(is_chsign(lp, rowno), value);
|
|
#ifdef DoMatrixRounding
|
|
value = roundToPrecision(value, mat->epsvalue);
|
|
if(value == 0) {
|
|
if((colnr > colnr1) || (rownr > rowno))
|
|
;
|
|
else
|
|
goto ForceDelete;
|
|
}
|
|
#endif
|
|
SET_MAT_ijA(jj_j, rowno, colnr1, value);
|
|
|
|
/* Adjust if we have inserted an element */
|
|
if((colnr > colnr1) || (rownr > rowno)) {
|
|
j--;
|
|
origidx--;
|
|
jj_j++;
|
|
delta++;
|
|
}
|
|
origidx++;
|
|
continue;
|
|
}
|
|
|
|
/* Shift the matrix element up by the active difference */
|
|
if(jj_j != origidx) {
|
|
COL_MAT_COPY(jj_j, origidx);
|
|
}
|
|
origidx++;
|
|
}
|
|
|
|
/* Update pending / incomplete column start position */
|
|
jj_j = origidx - j;
|
|
for(; k <= lendense; k++)
|
|
mat->col_end[k] = jj_j;
|
|
mat->row_end_valid = FALSE;
|
|
|
|
Done:
|
|
if(!isNZ)
|
|
FREE(row);
|
|
FREE(colno);
|
|
return( (MYBOOL) (newnz > 0) );
|
|
|
|
} /* mat_setrow */
|
|
#endif
|
|
|
|
STATIC int mat_appendrow(MATrec *mat, int count, REAL *row, int *colno, REAL mult, MYBOOL checkrowmode)
|
|
{
|
|
int i, j, jj = 0, stcol, elmnr, orignr, newnr, firstcol;
|
|
MYBOOL *addto = NULL, isA, isNZ;
|
|
REAL value, saved = 0;
|
|
lprec *lp = mat->lp;
|
|
|
|
/* Check if we are in row order mode and should add as column instead;
|
|
the matrix will be transposed at a later stage */
|
|
if(checkrowmode && mat->is_roworder)
|
|
return( mat_appendcol(mat, count, row, colno, mult, FALSE) );
|
|
|
|
/* Do initialization and validation */
|
|
isA = (MYBOOL) (mat == lp->matA);
|
|
isNZ = (MYBOOL) (colno != NULL);
|
|
if(isNZ && (count > 0)) {
|
|
if(count > 1)
|
|
sortREALByINT(row, colno, count, 0, TRUE);
|
|
if((colno[0] < 1) || (colno[count-1] > mat->columns))
|
|
return( 0 );
|
|
}
|
|
/* else if((row != NULL) && !mat->is_roworder) */
|
|
else if(!isNZ && (row != NULL) && !mat->is_roworder)
|
|
row[0] = 0;
|
|
|
|
/* Capture OF definition in row mode */
|
|
if(isA && mat->is_roworder) {
|
|
if(isNZ && (colno[0] == 0)) {
|
|
value = row[0];
|
|
#ifdef DoMatrixRounding
|
|
value = roundToPrecision(value, mat->epsvalue);
|
|
#endif
|
|
value = scaled_mat(lp, value, 0, lp->columns);
|
|
value = my_chsign(is_maxim(lp), value);
|
|
lp->orig_obj[lp->columns] = value;
|
|
count--;
|
|
row++;
|
|
colno++;
|
|
}
|
|
else if(!isNZ && (row != NULL) && (row[0] != 0)) {
|
|
value = saved = row[0];
|
|
#ifdef DoMatrixRounding
|
|
value = roundToPrecision(value, mat->epsvalue);
|
|
#endif
|
|
value = scaled_mat(lp, value, 0, lp->columns);
|
|
value = my_chsign(is_maxim(lp), value);
|
|
lp->orig_obj[lp->columns] = value;
|
|
row[0] = 0;
|
|
}
|
|
else
|
|
lp->orig_obj[lp->columns] = 0;
|
|
}
|
|
|
|
/* Optionally tally and map the new non-zero values */
|
|
firstcol = mat->columns + 1;
|
|
if(isNZ) {
|
|
newnr = count;
|
|
if(newnr) {
|
|
firstcol = colno[0];
|
|
jj = colno[newnr - 1];
|
|
}
|
|
}
|
|
else {
|
|
newnr = 0;
|
|
if(row != NULL) {
|
|
if(!allocMYBOOL(lp, &addto, mat->columns + 1, TRUE)) {
|
|
return( newnr );
|
|
}
|
|
for(i = mat->columns; i >= 1; i--) {
|
|
if(fabs(row[i]) > mat->epsvalue) {
|
|
addto[i] = TRUE;
|
|
firstcol = i;
|
|
newnr++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/* Make sure we have sufficient space */
|
|
if(!inc_mat_space(mat, newnr)) {
|
|
newnr = 0;
|
|
goto Done;
|
|
}
|
|
|
|
/* Insert the non-zero constraint values */
|
|
orignr = mat_nonzeros(mat) - 1;
|
|
elmnr = orignr + newnr;
|
|
|
|
for(j = mat->columns; j >= firstcol; j--) {
|
|
stcol = mat->col_end[j] - 1;
|
|
mat->col_end[j] = elmnr + 1;
|
|
|
|
/* Add a new non-zero entry */
|
|
if(((isNZ) && (j == jj)) || ((addto != NULL) && (addto[j]))) {
|
|
newnr--;
|
|
if(isNZ) {
|
|
value = row[newnr];
|
|
if(newnr)
|
|
jj = colno[newnr - 1];
|
|
else
|
|
jj = 0;
|
|
}
|
|
else
|
|
value = row[j];
|
|
#ifdef DoMatrixRounding
|
|
value = roundToPrecision(value, mat->epsvalue);
|
|
#endif
|
|
value *= mult;
|
|
if(isA) {
|
|
if(mat->is_roworder)
|
|
value = my_chsign(is_chsign(lp, j), value);
|
|
value = scaled_mat(lp, value, mat->rows, j);
|
|
}
|
|
SET_MAT_ijA(elmnr, mat->rows, j, value);
|
|
elmnr--;
|
|
}
|
|
|
|
/* Shift previous column entries down */
|
|
i = stcol - mat->col_end[j-1] + 1;
|
|
if(i > 0) {
|
|
orignr -= i;
|
|
elmnr -= i;
|
|
COL_MAT_MOVE(elmnr+1, orignr+1, i);
|
|
}
|
|
}
|
|
|
|
Done:
|
|
if(saved != 0)
|
|
row[0] = saved;
|
|
FREE(addto);
|
|
|
|
return( newnr );
|
|
|
|
}
|
|
|
|
STATIC int mat_appendcol(MATrec *mat, int count, REAL *column, int *rowno, REAL mult, MYBOOL checkrowmode)
|
|
{
|
|
int i, row, elmnr, lastnr;
|
|
REAL value;
|
|
MYBOOL isA, isNZ;
|
|
lprec *lp = mat->lp;
|
|
|
|
/* Check if we are in row order mode and should add as row instead;
|
|
the matrix will be transposed at a later stage */
|
|
if(checkrowmode && mat->is_roworder)
|
|
return( mat_appendrow(mat, count, column, rowno, mult, FALSE) );
|
|
|
|
/* Make sure we have enough space */
|
|
/*
|
|
if(!inc_mat_space(mat, mat->rows+1))
|
|
return( 0 );
|
|
*/
|
|
if(column == NULL)
|
|
i = 0;
|
|
else if(rowno != NULL)
|
|
i = count;
|
|
else {
|
|
int nrows = mat->rows;
|
|
|
|
elmnr = 0;
|
|
for(i = 1; i <= nrows; i++)
|
|
if(column[i] != 0)
|
|
elmnr++;
|
|
i = elmnr;
|
|
}
|
|
if((mat_nz_unused(mat) <= i) && !inc_mat_space(mat, i))
|
|
return( 0 );
|
|
|
|
/* Do initialization and validation */
|
|
isA = (MYBOOL) (mat == lp->matA);
|
|
isNZ = (MYBOOL) (column == NULL || rowno != NULL);
|
|
if(isNZ && (count > 0)) {
|
|
if(count > 1)
|
|
sortREALByINT(column, rowno, count, 0, TRUE);
|
|
if((rowno[0] < 0))
|
|
return( 0 );
|
|
}
|
|
if(rowno != NULL)
|
|
count--;
|
|
|
|
/* Append sparse regular constraint values */
|
|
elmnr = mat->col_end[mat->columns - 1];
|
|
if(column != NULL) {
|
|
row = -1;
|
|
for(i = ((isNZ || !mat->is_roworder) ? 0 : 1); i <= count ; i++) {
|
|
value = column[i];
|
|
if(fabs(value) > mat->epsvalue) {
|
|
if(isNZ) {
|
|
lastnr = row;
|
|
row = rowno[i];
|
|
/* Check if we have come to the Lagrangean constraints */
|
|
if(row > mat->rows)
|
|
break;
|
|
if(row <= lastnr)
|
|
return( -1 );
|
|
}
|
|
else
|
|
row = i;
|
|
#ifdef DoMatrixRounding
|
|
value = roundToPrecision(value, mat->epsvalue);
|
|
#endif
|
|
if(mat->is_roworder)
|
|
value *= mult;
|
|
else if(isA) {
|
|
value = my_chsign(is_chsign(lp, row), value);
|
|
value = scaled_mat(lp, value, row, mat->columns);
|
|
if(!mat->is_roworder && (row == 0)) {
|
|
lp->orig_obj[mat->columns] = value;
|
|
continue;
|
|
}
|
|
}
|
|
|
|
/* Store the item and update counters */
|
|
SET_MAT_ijA(elmnr, row, mat->columns, value);
|
|
elmnr++;
|
|
}
|
|
}
|
|
|
|
/* Fill dense Lagrangean constraints */
|
|
if(get_Lrows(lp) > 0)
|
|
mat_appendcol(lp->matL, get_Lrows(lp), column+mat->rows, NULL, mult, checkrowmode);
|
|
|
|
}
|
|
|
|
/* Set end of data */
|
|
mat->col_end[mat->columns] = elmnr;
|
|
|
|
return( mat->col_end[mat->columns] - mat->col_end[mat->columns-1] );
|
|
}
|
|
|
|
STATIC int mat_checkcounts(MATrec *mat, int *rownum, int *colnum, MYBOOL freeonexit)
|
|
{
|
|
int i, j, n;
|
|
int *rownr;
|
|
|
|
if(rownum == NULL)
|
|
allocINT(mat->lp, &rownum, mat->rows + 1, TRUE);
|
|
if(colnum == NULL)
|
|
allocINT(mat->lp, &colnum, mat->columns + 1, TRUE);
|
|
|
|
for(i = 1 ; i <= mat->columns; i++) {
|
|
j = mat->col_end[i - 1];
|
|
n = mat->col_end[i];
|
|
rownr = &COL_MAT_ROWNR(j);
|
|
for(; j < n;
|
|
j++, rownr += matRowColStep) {
|
|
colnum[i]++;
|
|
rownum[*rownr]++;
|
|
}
|
|
}
|
|
|
|
n = 0;
|
|
if((mat->lp->do_presolve != PRESOLVE_NONE) &&
|
|
(mat->lp->spx_trace || (mat->lp->verbose > NORMAL))) {
|
|
for(j = 1; j <= mat->columns; j++)
|
|
if(colnum[j] == 0) {
|
|
n++;
|
|
report(mat->lp, FULL, "mat_checkcounts: Variable %s is not used in any constraints\n",
|
|
get_col_name(mat->lp, j));
|
|
}
|
|
for(i = 0; i <= mat->rows; i++)
|
|
if(rownum[i] == 0) {
|
|
n++;
|
|
report(mat->lp, FULL, "mat_checkcounts: Constraint %s empty\n",
|
|
get_row_name(mat->lp, i));
|
|
}
|
|
}
|
|
|
|
if(freeonexit) {
|
|
FREE(rownum);
|
|
FREE(colnum);
|
|
}
|
|
|
|
return( n );
|
|
|
|
}
|
|
|
|
STATIC MYBOOL mat_validate(MATrec *mat)
|
|
/* Routine to make sure that row mapping arrays are valid */
|
|
{
|
|
int i, j, je, *rownum;
|
|
int *rownr, *colnr;
|
|
|
|
if(!mat->row_end_valid) {
|
|
|
|
MEMCLEAR(mat->row_end, mat->rows + 1);
|
|
allocINT(mat->lp, &rownum, mat->rows + 1, TRUE);
|
|
|
|
/* First tally row counts and then cumulate them */
|
|
j = mat_nonzeros(mat);
|
|
rownr = &COL_MAT_ROWNR(0);
|
|
for(i = 0; i < j; i++, rownr += matRowColStep)
|
|
mat->row_end[*rownr]++;
|
|
for(i = 1; i <= mat->rows; i++)
|
|
mat->row_end[i] += mat->row_end[i - 1];
|
|
|
|
/* Calculate the column index for every non-zero */
|
|
for(i = 1; i <= mat->columns; i++) {
|
|
j = mat->col_end[i - 1];
|
|
je = mat->col_end[i];
|
|
rownr = &COL_MAT_ROWNR(j);
|
|
colnr = &COL_MAT_COLNR(j);
|
|
for(; j < je; j++, rownr += matRowColStep, colnr += matRowColStep) {
|
|
#ifdef Paranoia
|
|
if(/*(*colnr < 0) || (*colnr > mat->columns) || (Normally violated in primal phase 1) */
|
|
(*rownr < 0) || (*rownr > mat->rows)) {
|
|
report(mat->lp, SEVERE, "mat_validate: Matrix value storage error row %d [0..%d], column %d [1..%d]\n",
|
|
*rownr, mat->rows, *colnr, mat->columns);
|
|
mat->lp->spx_status = UNKNOWNERROR;
|
|
return(FALSE);
|
|
}
|
|
#endif
|
|
*colnr = i;
|
|
if(*rownr == 0)
|
|
mat_set_rowmap(mat, rownum[*rownr],
|
|
*rownr, i, j);
|
|
else
|
|
mat_set_rowmap(mat, mat->row_end[*rownr - 1] + rownum[*rownr],
|
|
*rownr, i, j);
|
|
rownum[*rownr]++;
|
|
}
|
|
}
|
|
|
|
FREE(rownum);
|
|
mat->row_end_valid = TRUE;
|
|
}
|
|
|
|
if(mat == mat->lp->matA)
|
|
mat->lp->model_is_valid = TRUE;
|
|
return( TRUE );
|
|
}
|
|
|
|
MYBOOL mat_get_data(lprec *lp, int matindex, MYBOOL isrow, int **rownr, int **colnr, REAL **value)
|
|
{
|
|
MATrec *mat = lp->matA;
|
|
|
|
#if MatrixRowAccess == RAM_Index
|
|
if(isrow)
|
|
matindex = mat->row_mat[matindex];
|
|
if(rownr != NULL)
|
|
*rownr = &COL_MAT_ROWNR(matindex);
|
|
if(colnr != NULL)
|
|
*colnr = &COL_MAT_COLNR(matindex);
|
|
if(value != NULL)
|
|
*value = &COL_MAT_VALUE(matindex);
|
|
|
|
#else
|
|
if(isrow) {
|
|
if(rownr != NULL)
|
|
*rownr = &ROW_MAT_ROWNR(matindex);
|
|
if(colnr != NULL)
|
|
*colnr = &ROW_MAT_COLNR(matindex);
|
|
if(value != NULL)
|
|
*value = &ROW_MAT_VALUE(matindex);
|
|
}
|
|
else {
|
|
if(rownr != NULL)
|
|
*rownr = &COL_MAT_ROWNR(matindex);
|
|
if(colnr != NULL)
|
|
*colnr = &COL_MAT_COLNR(matindex);
|
|
if(value != NULL)
|
|
*value = &COL_MAT_VALUE(matindex);
|
|
}
|
|
|
|
#endif
|
|
|
|
return( TRUE );
|
|
}
|
|
|
|
|
|
MYBOOL mat_set_rowmap(MATrec *mat, int row_mat_index, int rownr, int colnr, int col_mat_index)
|
|
{
|
|
#if MatrixRowAccess == RAM_Index
|
|
mat->row_mat[row_mat_index] = col_mat_index;
|
|
|
|
#elif MatrixColAccess==CAM_Record
|
|
mat->row_mat[row_mat_index].rownr = rownr;
|
|
mat->row_mat[row_mat_index].colnr = colnr;
|
|
mat->row_mat[row_mat_index].value = COL_MAT_VALUE(col_mat_index);
|
|
|
|
#else /* if MatrixColAccess==CAM_Vector */
|
|
mat->row_mat_rownr[row_mat_index] = rownr;
|
|
mat->row_mat_colnr[row_mat_index] = colnr;
|
|
mat->row_mat_value[row_mat_index] = COL_MAT_VALUE(col_mat_index);
|
|
|
|
#endif
|
|
|
|
return( TRUE );
|
|
}
|
|
|
|
/* Implement combined binary/linear sub-search for matrix look-up */
|
|
int mat_findelm(MATrec *mat, int row, int column)
|
|
{
|
|
int low, high, mid, item;
|
|
|
|
#if 0
|
|
if(mat->row_end_valid && (row > 0) &&
|
|
(ROW_MAT_COLNR(mat->row_mat[(low = mat->row_end[row-1])]) == column))
|
|
return(low);
|
|
#endif
|
|
|
|
if((column < 1) || (column > mat->columns)) {
|
|
report(mat->lp, IMPORTANT, "mat_findelm: Column %d out of range\n", column);
|
|
return( -1 );
|
|
}
|
|
if((row < 0) || (row > mat->rows)) {
|
|
report(mat->lp, IMPORTANT, "mat_findelm: Row %d out of range\n", row);
|
|
return( -1 );
|
|
}
|
|
|
|
low = mat->col_end[column - 1];
|
|
high = mat->col_end[column] - 1;
|
|
if(low > high)
|
|
return( -2 );
|
|
|
|
/* Do binary search logic */
|
|
mid = (low+high) / 2;
|
|
item = COL_MAT_ROWNR(mid);
|
|
while(high - low > LINEARSEARCH) {
|
|
if(item < row) {
|
|
low = mid + 1;
|
|
mid = (low+high) / 2;
|
|
item = COL_MAT_ROWNR(mid);
|
|
}
|
|
else if(item > row) {
|
|
high = mid - 1;
|
|
mid = (low+high) / 2;
|
|
item = COL_MAT_ROWNR(mid);
|
|
}
|
|
else {
|
|
low = mid;
|
|
high = mid;
|
|
}
|
|
}
|
|
|
|
/* Do linear scan search logic */
|
|
if((high > low) && (high - low <= LINEARSEARCH)) {
|
|
item = COL_MAT_ROWNR(low);
|
|
while((low < high) && (item < row)) {
|
|
low++;
|
|
item = COL_MAT_ROWNR(low);
|
|
}
|
|
if(item == row)
|
|
high = low;
|
|
}
|
|
|
|
if((low == high) && (row == item))
|
|
return( low );
|
|
else
|
|
return( -2 );
|
|
}
|
|
|
|
int mat_findins(MATrec *mat, int row, int column, int *insertpos, MYBOOL validate)
|
|
{
|
|
int low, high, mid, item, exitvalue, insvalue;
|
|
|
|
#if 0
|
|
if(mat->row_end_valid && (row > 0) &&
|
|
(ROW_MAT_COLNR(mat->row_mat[(low = mat->row_end[row-1])]) == column)) {
|
|
insvalue = low;
|
|
exitvalue = low;
|
|
goto Done;
|
|
}
|
|
#endif
|
|
|
|
insvalue = -1;
|
|
|
|
if((column < 1) || (column > mat->columns)) {
|
|
if((column > 0) && !validate) {
|
|
insvalue = mat->col_end[mat->columns];
|
|
exitvalue = -2;
|
|
goto Done;
|
|
}
|
|
report(mat->lp, IMPORTANT, "mat_findins: Column %d out of range\n", column);
|
|
exitvalue = -1;
|
|
goto Done;
|
|
}
|
|
if((row < 0) || (row > mat->rows)) {
|
|
if((row >= 0) && !validate) {
|
|
insvalue = mat->col_end[column];
|
|
exitvalue = -2;
|
|
goto Done;
|
|
}
|
|
report(mat->lp, IMPORTANT, "mat_findins: Row %d out of range\n", row);
|
|
exitvalue = -1;
|
|
goto Done;
|
|
}
|
|
|
|
low = mat->col_end[column - 1];
|
|
insvalue = low;
|
|
high = mat->col_end[column] - 1;
|
|
if(low > high) {
|
|
exitvalue = -2;
|
|
goto Done;
|
|
}
|
|
|
|
/* Do binary search logic */
|
|
mid = (low+high) / 2;
|
|
item = COL_MAT_ROWNR(mid);
|
|
while(high - low > LINEARSEARCH) {
|
|
if(item < row) {
|
|
low = mid + 1;
|
|
mid = (low+high) / 2;
|
|
item = COL_MAT_ROWNR(mid);
|
|
}
|
|
else if(item > row) {
|
|
high = mid - 1;
|
|
mid = (low+high) / 2;
|
|
item = COL_MAT_ROWNR(mid);
|
|
}
|
|
else {
|
|
low = mid;
|
|
high = mid;
|
|
}
|
|
}
|
|
|
|
/* Do linear scan search logic */
|
|
if((high > low) && (high - low <= LINEARSEARCH)) {
|
|
item = COL_MAT_ROWNR(low);
|
|
while((low < high) && (item < row)) {
|
|
low++;
|
|
item = COL_MAT_ROWNR(low);
|
|
}
|
|
if(item == row)
|
|
high = low;
|
|
}
|
|
|
|
insvalue = low;
|
|
if((low == high) && (row == item))
|
|
exitvalue = low;
|
|
else {
|
|
if((low < mat->col_end[column]) && (COL_MAT_ROWNR(low) < row))
|
|
insvalue++;
|
|
exitvalue = -2;
|
|
}
|
|
|
|
Done:
|
|
if(insertpos != NULL)
|
|
(*insertpos) = insvalue;
|
|
return( exitvalue );
|
|
}
|
|
|
|
STATIC REAL mat_getitem(MATrec *mat, int row, int column)
|
|
{
|
|
int elmnr;
|
|
|
|
#ifdef DirectOverrideOF
|
|
if((row == 0) && (mat == mat->lp->matA) && (mat->lp->OF_override != NULL))
|
|
return( mat->lp->OF_override[column] );
|
|
else
|
|
#endif
|
|
{
|
|
elmnr = mat_findelm(mat, row, column);
|
|
if(elmnr >= 0)
|
|
return( COL_MAT_VALUE(elmnr) );
|
|
else
|
|
return( 0 );
|
|
}
|
|
}
|
|
|
|
STATIC MYBOOL mat_additem(MATrec *mat, int row, int column, REAL delta)
|
|
{
|
|
int elmnr;
|
|
|
|
#ifdef DirectOverrideOF
|
|
if((row == 0) && (mat == mat->lp->matA) && (mat->lp->OF_override != NULL))
|
|
return( mat->lp->OF_override[column] );
|
|
else
|
|
#endif
|
|
{
|
|
elmnr = mat_findelm(mat, row, column);
|
|
if(elmnr >= 0) {
|
|
COL_MAT_VALUE(elmnr) += delta;
|
|
return( TRUE );
|
|
}
|
|
else {
|
|
mat_setitem(mat, row, column, delta);
|
|
return( FALSE );
|
|
}
|
|
}
|
|
}
|
|
|
|
STATIC MYBOOL mat_setitem(MATrec *mat, int row, int column, REAL value)
|
|
{
|
|
return( mat_setvalue(mat, row, column, value, FALSE) );
|
|
}
|
|
|
|
STATIC void mat_multrow(MATrec *mat, int row_nr, REAL mult)
|
|
{
|
|
int i, k1, k2;
|
|
|
|
#if 0
|
|
if(row_nr == 0) {
|
|
k2 = mat->col_end[0];
|
|
for(i = 1; i <= mat->columns; i++) {
|
|
k1 = k2;
|
|
k2 = mat->col_end[i];
|
|
if((k1 < k2) && (COL_MAT_ROWNR(k1) == row_nr))
|
|
COL_MAT_VALUE(k1) *= mult;
|
|
}
|
|
}
|
|
else if(mat_validate(mat)) {
|
|
if(row_nr == 0)
|
|
k1 = 0;
|
|
else
|
|
#else
|
|
if(mat_validate(mat)) {
|
|
if(row_nr == 0)
|
|
k1 = 0;
|
|
else
|
|
#endif
|
|
k1 = mat->row_end[row_nr-1];
|
|
k2 = mat->row_end[row_nr];
|
|
for(i = k1; i < k2; i++)
|
|
ROW_MAT_VALUE(i) *= mult;
|
|
}
|
|
}
|
|
|
|
STATIC void mat_multcol(MATrec *mat, int col_nr, REAL mult, MYBOOL DoObj)
|
|
{
|
|
int i, ie;
|
|
MYBOOL isA;
|
|
|
|
#ifdef Paranoia
|
|
if((col_nr < 1) || (col_nr > mat->columns)) {
|
|
report(mat->lp, IMPORTANT, "mult_column: Column %d out of range\n", col_nr);
|
|
return;
|
|
}
|
|
#endif
|
|
if(mult == 1.0)
|
|
return;
|
|
|
|
isA = (MYBOOL) (mat == mat->lp->matA);
|
|
|
|
ie = mat->col_end[col_nr];
|
|
for(i = mat->col_end[col_nr - 1]; i < ie; i++)
|
|
COL_MAT_VALUE(i) *= mult;
|
|
if(isA) {
|
|
if(DoObj)
|
|
mat->lp->orig_obj[col_nr] *= mult;
|
|
if(get_Lrows(mat->lp) > 0)
|
|
mat_multcol(mat->lp->matL, col_nr, mult, DoObj);
|
|
}
|
|
}
|
|
|
|
STATIC void mat_multadd(MATrec *mat, REAL *lhsvector, int varnr, REAL mult)
|
|
{
|
|
int colnr;
|
|
register int ib, ie, *matRownr;
|
|
register REAL *matValue;
|
|
|
|
/* Handle case of a slack variable */
|
|
if(varnr <= mat->lp->rows) {
|
|
lhsvector[varnr] += mult;
|
|
return;
|
|
}
|
|
|
|
/* Do operation on the objective */
|
|
if(mat->lp->matA == mat)
|
|
lhsvector[0] += get_OF_active(mat->lp, varnr, mult);
|
|
|
|
/* Scan the constraint matrix target columns */
|
|
colnr = varnr - mat->lp->rows;
|
|
ib = mat->col_end[colnr - 1];
|
|
ie = mat->col_end[colnr];
|
|
if(ib < ie) {
|
|
|
|
/* Initialize pointers */
|
|
matRownr = &COL_MAT_ROWNR(ib);
|
|
matValue = &COL_MAT_VALUE(ib);
|
|
|
|
/* Then loop over all regular rows */
|
|
for(; ib < ie;
|
|
ib++, matValue += matValueStep, matRownr += matRowColStep) {
|
|
lhsvector[*matRownr] += mult * (*matValue);
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
STATIC MYBOOL mat_setvalue(MATrec *mat, int Row, int Column, REAL Value, MYBOOL doscale)
|
|
{
|
|
int elmnr, lastelm, i, RowA = Row, ColumnA = Column;
|
|
MYBOOL isA;
|
|
|
|
/* This function is inefficient if used to add new matrix entries in
|
|
other places than at the end of the matrix. OK for replacing existing
|
|
a non-zero value with another non-zero value */
|
|
isA = (MYBOOL) (mat == mat->lp->matA);
|
|
if(mat->is_roworder)
|
|
swapINT(&Row, &Column);
|
|
|
|
/* Set small numbers to zero */
|
|
if(fabs(Value) < mat->epsvalue)
|
|
Value = 0;
|
|
#ifdef DoMatrixRounding
|
|
else
|
|
Value = roundToPrecision(Value, mat->epsvalue);
|
|
#endif
|
|
|
|
/* Check if we need to update column space */
|
|
if(Column > mat->columns) {
|
|
if(isA)
|
|
inc_col_space(mat->lp, ColumnA - mat->columns);
|
|
else
|
|
inc_matcol_space(mat, Column - mat->columns);
|
|
}
|
|
|
|
/* Find out if we already have such an entry, or return insertion point */
|
|
i = mat_findins(mat, Row, Column, &elmnr, FALSE);
|
|
if(i == -1)
|
|
return(FALSE);
|
|
|
|
if(isA)
|
|
set_action(&mat->lp->spx_action, ACTION_REBASE | ACTION_RECOMPUTE | ACTION_REINVERT);
|
|
|
|
if(i >= 0) {
|
|
/* there is an existing entry */
|
|
if(fabs(Value) > mat->epsvalue) { /* we replace it by something non-zero */
|
|
if(isA) {
|
|
Value = my_chsign(is_chsign(mat->lp, RowA), Value);
|
|
if(doscale && mat->lp->scaling_used)
|
|
Value = scaled_mat(mat->lp, Value, RowA, ColumnA);
|
|
}
|
|
COL_MAT_VALUE(elmnr) = Value;
|
|
}
|
|
else { /* setting existing non-zero entry to zero. Remove the entry */
|
|
/* This might remove an entire column, or leave just a bound. No
|
|
nice solution for that yet */
|
|
|
|
/* Shift up tail end of the matrix */
|
|
lastelm = mat_nonzeros(mat);
|
|
#if 0
|
|
for(i = elmnr; i < lastelm ; i++) {
|
|
COL_MAT_COPY(i, i + 1);
|
|
}
|
|
#else
|
|
lastelm -= elmnr;
|
|
COL_MAT_MOVE(elmnr, elmnr + 1, lastelm);
|
|
#endif
|
|
for(i = Column; i <= mat->columns; i++)
|
|
mat->col_end[i]--;
|
|
|
|
mat->row_end_valid = FALSE;
|
|
}
|
|
}
|
|
else if(fabs(Value) > mat->epsvalue) {
|
|
/* no existing entry. make new one only if not nearly zero */
|
|
/* check if more space is needed for matrix */
|
|
if(!inc_mat_space(mat, 1))
|
|
return(FALSE);
|
|
|
|
if(Column > mat->columns) {
|
|
i = mat->columns + 1;
|
|
if(isA)
|
|
shift_coldata(mat->lp, i, ColumnA - mat->columns, NULL);
|
|
else
|
|
mat_shiftcols(mat, &i, Column - mat->columns, NULL);
|
|
}
|
|
|
|
/* Shift down tail end of the matrix by one */
|
|
lastelm = mat_nonzeros(mat);
|
|
#if 1 /* Does compiler optimization work better here? */
|
|
for(i = lastelm; i > elmnr ; i--) {
|
|
COL_MAT_COPY(i, i - 1);
|
|
}
|
|
#else
|
|
lastelm -= elmnr - 1;
|
|
COL_MAT_MOVE(elmnr + 1, elmnr, lastelm);
|
|
#endif
|
|
|
|
/* Set new element */
|
|
if(isA) {
|
|
Value = my_chsign(is_chsign(mat->lp, RowA), Value);
|
|
if(doscale)
|
|
Value = scaled_mat(mat->lp, Value, RowA, ColumnA);
|
|
}
|
|
SET_MAT_ijA(elmnr, Row, Column, Value);
|
|
|
|
/* Update column indexes */
|
|
for(i = Column; i <= mat->columns; i++)
|
|
mat->col_end[i]++;
|
|
|
|
mat->row_end_valid = FALSE;
|
|
}
|
|
|
|
if(isA && (mat->lp->var_is_free != NULL) && (mat->lp->var_is_free[ColumnA] > 0))
|
|
return( mat_setvalue(mat, RowA, mat->lp->var_is_free[ColumnA], -Value, doscale) );
|
|
return(TRUE);
|
|
}
|
|
|
|
STATIC MYBOOL mat_appendvalue(MATrec *mat, int Row, REAL Value)
|
|
{
|
|
int *elmnr, Column = mat->columns;
|
|
|
|
/* Set small numbers to zero */
|
|
if(fabs(Value) < mat->epsvalue)
|
|
Value = 0;
|
|
#ifdef DoMatrixRounding
|
|
else
|
|
Value = roundToPrecision(Value, mat->epsvalue);
|
|
#endif
|
|
|
|
/* Check if more space is needed for matrix */
|
|
if(!inc_mat_space(mat, 1))
|
|
return(FALSE);
|
|
|
|
#ifdef Paranoia
|
|
/* Check valid indeces */
|
|
if((Row < 0) || (Row > mat->rows)) {
|
|
report(mat->lp, SEVERE, "mat_appendvalue: Invalid row index %d specified\n", Row);
|
|
return(FALSE);
|
|
}
|
|
#endif
|
|
|
|
/* Get insertion point and set value */
|
|
elmnr = mat->col_end + Column;
|
|
SET_MAT_ijA((*elmnr), Row, Column, Value);
|
|
|
|
/* Update column count */
|
|
(*elmnr)++;
|
|
mat->row_end_valid = FALSE;
|
|
|
|
return(TRUE);
|
|
}
|
|
|
|
STATIC MYBOOL mat_equalRows(MATrec *mat, int baserow, int comprow)
|
|
{
|
|
MYBOOL status = FALSE;
|
|
|
|
if(mat_validate(mat)) {
|
|
int bj1 = 0, ej1, bj2 = 0, ej2;
|
|
|
|
/* Get starting and ending positions */
|
|
if(baserow >= 0)
|
|
bj1 = mat->row_end[baserow-1];
|
|
ej1 = mat->row_end[baserow];
|
|
if(comprow >= 0)
|
|
bj2 = mat->row_end[comprow-1];
|
|
ej2 = mat->row_end[comprow];
|
|
/* Fail if row lengths are unequal */
|
|
if((ej1-bj1) != (ej2-bj2))
|
|
return( status );
|
|
|
|
/* Compare column index and value, element by element */
|
|
for(; bj1 < ej1; bj1++, bj2++) {
|
|
if(COL_MAT_COLNR(bj1) != COL_MAT_COLNR(bj2))
|
|
break;
|
|
#if 1
|
|
if(fabs(get_mat_byindex(mat->lp, bj1, TRUE, FALSE)-get_mat_byindex(mat->lp, bj2, TRUE, FALSE)) > mat->lp->epsprimal)
|
|
#else
|
|
if(fabs(COL_MAT_VALUE(bj1)-COL_MAT_VALUE(bj2)) > mat->lp->epsprimal)
|
|
#endif
|
|
break;
|
|
}
|
|
status = (MYBOOL) (bj1 == ej1);
|
|
}
|
|
return( status );
|
|
}
|
|
|
|
STATIC int mat_findcolumn(MATrec *mat, int matindex)
|
|
{
|
|
int j;
|
|
|
|
for(j = 1; j <= mat->columns; j++) {
|
|
if(matindex < mat->col_end[j])
|
|
break;
|
|
}
|
|
return(j);
|
|
}
|
|
|
|
STATIC int mat_expandcolumn(MATrec *mat, int colnr, REAL *column, int *nzlist, MYBOOL signedA)
|
|
{
|
|
MYBOOL isA = (MYBOOL) (mat->lp->matA == mat);
|
|
int i, ie, j, nzcount = 0;
|
|
REAL *matValue;
|
|
int *matRownr;
|
|
|
|
signedA &= isA;
|
|
|
|
/* Retrieve a column from the user data matrix A */
|
|
MEMCLEAR(column, mat->rows + 1);
|
|
if(isA) {
|
|
column[0] = mat->lp->orig_obj[colnr];
|
|
if(signedA && is_chsign(mat->lp, 0))
|
|
column[0] = -column[0];
|
|
}
|
|
|
|
i = mat->col_end[colnr - 1];
|
|
ie = mat->col_end[colnr];
|
|
matRownr = &COL_MAT_ROWNR(i);
|
|
matValue = &COL_MAT_VALUE(i);
|
|
for(; i < ie;
|
|
i++, matRownr += matRowColStep, matValue += matValueStep) {
|
|
j = *matRownr;
|
|
column[j] = *matValue;
|
|
if(signedA && is_chsign(mat->lp, j))
|
|
column[j] = -column[j];
|
|
nzcount++;
|
|
if(nzlist != NULL)
|
|
nzlist[nzcount] = j;
|
|
}
|
|
if(nzlist != NULL)
|
|
nzlist[0] = nzcount;
|
|
return( nzcount );
|
|
}
|
|
|
|
STATIC MYBOOL mat_computemax(MATrec *mat)
|
|
{
|
|
int *rownr = &COL_MAT_ROWNR(0),
|
|
*colnr = &COL_MAT_COLNR(0),
|
|
i = 0, ie = mat->col_end[mat->columns], ez = 0;
|
|
REAL *value = &COL_MAT_VALUE(0), epsmachine = mat->lp->epsmachine, absvalue;
|
|
|
|
/* Prepare arrays */
|
|
if(!allocREAL(mat->lp, &mat->colmax, mat->columns_alloc+1, AUTOMATIC) ||
|
|
!allocREAL(mat->lp, &mat->rowmax, mat->rows_alloc+1, AUTOMATIC))
|
|
return( FALSE );
|
|
MEMCLEAR(mat->colmax, mat->columns+1);
|
|
MEMCLEAR(mat->rowmax, mat->rows+1);
|
|
|
|
/* Obtain the row and column maxima in one sweep */
|
|
mat->dynrange = mat->lp->infinite;
|
|
for(; i < ie;
|
|
i++, rownr += matRowColStep, colnr += matRowColStep, value += matValueStep) {
|
|
absvalue = fabs(*value);
|
|
SETMAX(mat->colmax[*colnr], absvalue);
|
|
SETMAX(mat->rowmax[*rownr], absvalue);
|
|
SETMIN(mat->dynrange, absvalue);
|
|
if(absvalue < epsmachine)
|
|
ez++;
|
|
}
|
|
|
|
/* Lastly, compute the global maximum and get the dynamic range */
|
|
for(i = 1; i <= mat->rows; i++)
|
|
SETMAX(mat->rowmax[0], mat->rowmax[i]);
|
|
mat->infnorm = mat->colmax[0] = mat->rowmax[0];
|
|
if(mat->dynrange == 0) {
|
|
report(mat->lp, SEVERE, "%d matrix contains zero-valued coefficients.\n", ez);
|
|
mat->dynrange = mat->lp->infinite;
|
|
}
|
|
else {
|
|
mat->dynrange = mat->infnorm / mat->dynrange;
|
|
if(ez > 0)
|
|
report(mat->lp, IMPORTANT, "%d matrix coefficients below machine precision were found.\n", ez);
|
|
}
|
|
|
|
return( TRUE );
|
|
}
|
|
|
|
STATIC MYBOOL mat_transpose(MATrec *mat)
|
|
{
|
|
int i, j, nz, k;
|
|
MYBOOL status;
|
|
|
|
status = mat_validate(mat);
|
|
if(status) {
|
|
|
|
/* Create a column-ordered sparse element list; "column" index must be shifted */
|
|
nz = mat_nonzeros(mat);
|
|
if(nz > 0) {
|
|
#if MatrixColAccess==CAM_Record
|
|
MATitem *newmat;
|
|
newmat = (MATitem *) malloc((mat->mat_alloc) * sizeof(*(mat->col_mat)));
|
|
j = mat->row_end[0];
|
|
for(i = nz-1; i >= j ; i--) {
|
|
k = i-j;
|
|
newmat[k] = mat->col_mat[mat->row_mat[i]];
|
|
newmat[k].row_nr = newmat[k].col_nr;
|
|
}
|
|
for(i = j-1; i >= 0 ; i--) {
|
|
k = nz-j+i;
|
|
newmat[k] = mat->col_mat[mat->row_mat[i]];
|
|
newmat[k].row_nr = newmat[k].col_nr;
|
|
}
|
|
swapPTR((void **) &mat->col_mat, (void **) &newmat);
|
|
FREE(newmat);
|
|
#else /*if MatrixColAccess==CAM_Vector*/
|
|
REAL *newValue = NULL;
|
|
int *newRownr = NULL;
|
|
allocREAL(mat->lp, &newValue, mat->mat_alloc, FALSE);
|
|
allocINT(mat->lp, &newRownr, mat->mat_alloc, FALSE);
|
|
|
|
j = mat->row_end[0];
|
|
for(i = nz-1; i >= j ; i--) {
|
|
k = i-j;
|
|
newValue[k] = ROW_MAT_VALUE(i);
|
|
newRownr[k] = ROW_MAT_COLNR(i);
|
|
}
|
|
for(i = j-1; i >= 0 ; i--) {
|
|
k = nz-j+i;
|
|
newValue[k] = ROW_MAT_VALUE(i);
|
|
newRownr[k] = ROW_MAT_COLNR(i);
|
|
}
|
|
|
|
swapPTR((void **) &mat->col_mat_rownr, (void **) &newRownr);
|
|
swapPTR((void **) &mat->col_mat_value, (void **) &newValue);
|
|
FREE(newValue);
|
|
FREE(newRownr);
|
|
#endif
|
|
}
|
|
|
|
/* Transfer row start to column start position; must adjust for different offsets */
|
|
if(mat->rows == mat->rows_alloc)
|
|
inc_matcol_space(mat, 1);
|
|
j = mat->row_end[0];
|
|
for(i = mat->rows; i >= 1; i--)
|
|
mat->row_end[i] -= j;
|
|
mat->row_end[mat->rows] = nz;
|
|
swapPTR((void **) &mat->row_end, (void **) &mat->col_end);
|
|
|
|
/* Swap arrays of maximum values */
|
|
swapPTR((void **) &mat->rowmax, (void **) &mat->colmax);
|
|
|
|
/* Swap array sizes */
|
|
swapINT(&mat->rows, &mat->columns);
|
|
swapINT(&mat->rows_alloc, &mat->columns_alloc);
|
|
|
|
/* Finally set current storage mode */
|
|
mat->is_roworder = (MYBOOL) !mat->is_roworder;
|
|
mat->row_end_valid = FALSE;
|
|
}
|
|
return(status);
|
|
}
|
|
|
|
|
|
/* ---------------------------------------------------------------------------------- */
|
|
/* Change-tracking routines */
|
|
/* ---------------------------------------------------------------------------------- */
|
|
STATIC DeltaVrec *createUndoLadder(lprec *lp, int levelitems, int maxlevels)
|
|
{
|
|
DeltaVrec *hold;
|
|
|
|
hold = (DeltaVrec *) malloc(sizeof(*hold));
|
|
hold->lp = lp;
|
|
hold->activelevel = 0;
|
|
hold->tracker = mat_create(lp, levelitems, 0, 0.0);
|
|
inc_matcol_space(hold->tracker, maxlevels);
|
|
return( hold );
|
|
}
|
|
STATIC int incrementUndoLadder(DeltaVrec *DV)
|
|
{
|
|
DV->activelevel++;
|
|
inc_matcol_space(DV->tracker, 1);
|
|
mat_shiftcols(DV->tracker, &(DV->activelevel), 1, NULL);
|
|
DV->tracker->columns++;
|
|
return(DV->activelevel);
|
|
}
|
|
STATIC MYBOOL modifyUndoLadder(DeltaVrec *DV, int itemno, REAL target[], REAL newvalue)
|
|
{
|
|
MYBOOL status;
|
|
int varindex = itemno;
|
|
REAL oldvalue = target[itemno];
|
|
|
|
#ifndef UseMilpSlacksRCF /* Check if we should include ranged constraints */
|
|
varindex -= DV->lp->rows;
|
|
#endif
|
|
status = mat_appendvalue(DV->tracker, varindex, oldvalue);
|
|
target[itemno] = newvalue;
|
|
return(status);
|
|
}
|
|
STATIC int countsUndoLadder(DeltaVrec *DV)
|
|
{
|
|
if(DV->activelevel > 0)
|
|
return( mat_collength(DV->tracker, DV->activelevel) );
|
|
else
|
|
return( 0 );
|
|
}
|
|
STATIC int restoreUndoLadder(DeltaVrec *DV, REAL target[])
|
|
{
|
|
int iD = 0;
|
|
|
|
if(DV->activelevel > 0) {
|
|
MATrec *mat = DV->tracker;
|
|
int iB = mat->col_end[DV->activelevel-1],
|
|
iE = mat->col_end[DV->activelevel],
|
|
*matRownr = &COL_MAT_ROWNR(iB);
|
|
REAL *matValue = &COL_MAT_VALUE(iB),
|
|
oldvalue;
|
|
|
|
/* Restore the values */
|
|
iD = iE-iB;
|
|
for(; iB < iE; iB++, matValue += matValueStep, matRownr += matRowColStep) {
|
|
oldvalue = *matValue;
|
|
#ifdef UseMilpSlacksRCF /* Check if we should include ranged constraints */
|
|
target[(*matRownr)] = oldvalue;
|
|
#else
|
|
target[DV->lp->rows+(*matRownr)] = oldvalue;
|
|
#endif
|
|
}
|
|
|
|
/* Get rid of the changes */
|
|
mat_shiftcols(DV->tracker, &(DV->activelevel), -1, NULL);
|
|
}
|
|
|
|
return(iD);
|
|
}
|
|
STATIC int decrementUndoLadder(DeltaVrec *DV)
|
|
{
|
|
int deleted = 0;
|
|
|
|
if(DV->activelevel > 0) {
|
|
deleted = mat_shiftcols(DV->tracker, &(DV->activelevel), -1, NULL);
|
|
DV->activelevel--;
|
|
DV->tracker->columns--;
|
|
}
|
|
return(deleted);
|
|
}
|
|
STATIC MYBOOL freeUndoLadder(DeltaVrec **DV)
|
|
{
|
|
if((DV == NULL) || (*DV == NULL))
|
|
return(FALSE);
|
|
|
|
mat_free(&((*DV)->tracker));
|
|
FREE(*DV);
|
|
return(TRUE);
|
|
}
|
|
|
|
STATIC MYBOOL appendUndoPresolve(lprec *lp, MYBOOL isprimal, REAL beta, int colnrDep)
|
|
{
|
|
MATrec *mat;
|
|
|
|
/* Point to correct undo structure */
|
|
if(isprimal)
|
|
mat = lp->presolve_undo->primalundo->tracker;
|
|
else
|
|
mat = lp->presolve_undo->dualundo->tracker;
|
|
|
|
/* Append the data */
|
|
if((colnrDep > 0) && (beta != 0) &&
|
|
(mat != NULL) && (mat->col_tag[0] > 0)) {
|
|
int ix = mat->col_tag[0];
|
|
#if 0
|
|
report(lp, NORMAL, "appendUndoPresolve: %s %g * x%d\n",
|
|
( beta < 0 ? "-" : "+"), fabs(beta), colnrDep);
|
|
#endif
|
|
|
|
/* Do normal user variable case */
|
|
if(colnrDep <= lp->columns)
|
|
mat_setvalue(mat, colnrDep, ix, beta, FALSE);
|
|
|
|
/* Handle case where a slack variable is referenced */
|
|
else {
|
|
int ipos, jx = mat->col_tag[ix];
|
|
mat_setvalue(mat, jx, ix, beta, FALSE);
|
|
jx = mat_findins(mat, jx, ix, &ipos, FALSE);
|
|
COL_MAT_ROWNR(ipos) = colnrDep;
|
|
}
|
|
return( TRUE );
|
|
}
|
|
else
|
|
return( FALSE );
|
|
}
|
|
STATIC MYBOOL addUndoPresolve(lprec *lp, MYBOOL isprimal, int colnrElim, REAL alpha, REAL beta, int colnrDep)
|
|
{
|
|
int ix;
|
|
DeltaVrec **DV;
|
|
MATrec *mat;
|
|
presolveundorec *psdata = lp->presolve_undo;
|
|
|
|
/* Point to and initialize undo structure at first call */
|
|
if(isprimal) {
|
|
DV = &(psdata->primalundo);
|
|
if(*DV == NULL) {
|
|
*DV = createUndoLadder(lp, lp->columns+1, lp->columns);
|
|
mat = (*DV)->tracker;
|
|
mat->epsvalue = lp->matA->epsvalue;
|
|
allocINT(lp, &(mat->col_tag), lp->columns+1, FALSE);
|
|
mat->col_tag[0] = 0;
|
|
}
|
|
}
|
|
else {
|
|
DV = &(psdata->dualundo);
|
|
if(*DV == NULL) {
|
|
*DV = createUndoLadder(lp, lp->rows+1, lp->rows);
|
|
mat = (*DV)->tracker;
|
|
mat->epsvalue = lp->matA->epsvalue;
|
|
allocINT(lp, &(mat->col_tag), lp->rows+1, FALSE);
|
|
mat->col_tag[0] = 0;
|
|
}
|
|
}
|
|
mat = (*DV)->tracker;
|
|
#if 0
|
|
report(lp, NORMAL, "addUndoPresolve: x%d = %g %s %g * x%d\n",
|
|
colnrElim, alpha, ( beta < 0 ? "-" : "+"), fabs(beta), colnrDep);
|
|
#endif
|
|
/* Add the data */
|
|
ix = mat->col_tag[0] = incrementUndoLadder(*DV);
|
|
mat->col_tag[ix] = colnrElim;
|
|
if(alpha != 0)
|
|
mat_setvalue(mat, 0, ix, alpha, FALSE);
|
|
/* mat_appendvalue(*mat, 0, alpha);*/
|
|
if((colnrDep > 0) && (beta != 0)) {
|
|
if(colnrDep > lp->columns)
|
|
return( appendUndoPresolve(lp, isprimal, beta, colnrDep) );
|
|
else
|
|
mat_setvalue(mat, colnrDep, ix, beta, FALSE);
|
|
}
|
|
|
|
return( TRUE );
|
|
}
|
|
|
|
|
|
|
|
/* ---------------------------------------------------------------------------------- */
|
|
/* High level matrix inverse and product routines in lp_solve */
|
|
/* ---------------------------------------------------------------------------------- */
|
|
|
|
/* ---------------------------------------------------------------------------------- */
|
|
/* A brief description of the basis inverse and factorization logic in lp_solve */
|
|
/* ---------------------------------------------------------------------------------- */
|
|
/*
|
|
|
|
In order to better understand the legacy code for operating with the
|
|
basis and its factorization in lp_solve I (KE) will briefly explain
|
|
the conventions and associated matrix algebra. Note that with lp_solve
|
|
version 5.5, it is also possible to direct lp_solve to use the traditional
|
|
(textbook) format by setting the obj_in_B parameter to FALSE.
|
|
|
|
The matrix description of a linear program (as represented by lp_solve) goes
|
|
like this:
|
|
|
|
maximize c'x
|
|
subject to r <= Ax <= b
|
|
where l <= x <= u
|
|
|
|
The matrix A is partitioned into two column sets [B|N], where B is
|
|
a square matrix of "basis" variables containing non-fixed
|
|
variables of the linear program at any given stage and N is the
|
|
submatrix of corresponding non-basic, fixed variables. The
|
|
variables (columns) in N may be fixed at their lower or upper levels.
|
|
|
|
Similarly, the c vector is partitioned into the basic and non-basic
|
|
parts [z|n].
|
|
|
|
While lp_solve stores the objective vector c in a dense format, and
|
|
the constraint matrix A in a (fairly standard) sparse format, the
|
|
column vectors passed to the factorization routine include the
|
|
objective coefficient at row index 0. (In versions of lp_solve
|
|
before v5.2, c was actually explicitly stored as the 0-th row of A).
|
|
The expanded matrix may be called the "A~" form and looks like this:
|
|
|
|
A~ = [ c ]
|
|
[ A ]
|
|
|
|
Linear programming involves solving linear equations based on the
|
|
square basis matrix B, which includes is a subset of columns from A~.
|
|
The implications of the common storage of c and A (e.g. A~) vs. the
|
|
inverse / factorization of B for the operations and updates performed
|
|
by the simplex routine therefore needs to be understood. As a consquence
|
|
of A~, in lp_solve B is stored in an expanded, bordered format using the
|
|
following (non-singular) representation:
|
|
|
|
B~ = [ 1 z ]
|
|
[ 0 B ]
|
|
|
|
Any basis inversion / factorization engine used by lp_solve must therefore
|
|
explicitly represent and handle the implications of this structure for
|
|
associated matrix operations.
|
|
|
|
The standard matrix formula for computing the inverse of a bordered
|
|
matrix shows what the inversion of B~ actually produces:
|
|
|
|
Inv(B~) = [ 1 -z*Inv(B) ]
|
|
[ 0 Inv(B) ]
|
|
|
|
The A~ and B~ representations require awareness by the developer of the side
|
|
effects of the presence of the top row when doing product operations such as
|
|
b'N, btran and ftran. Note in particular z*Inv(B) in the top row of Inv(B~),
|
|
which is actually the dual solution vector of the given basis. This fact
|
|
makes a very common update in the simplex algorithm (reduced costs) returnable
|
|
as a vector simply by setting 1 at the top of a vector being pre-multiplied
|
|
with Inv(B~).
|
|
|
|
However, if the objective vector (c) is changed, the expanded representation
|
|
requires that B / B~ be refactorized. Also, when doing FTRAN, BTRAN
|
|
and x'A-type operations, you will patently get the incorrect result
|
|
if you simply copy the operations given in textbooks. First I'll show the
|
|
results of an FTRAN operation:
|
|
|
|
Bx = a ==> x = FTRAN(a)
|
|
|
|
In lp_solve, this operation solves:
|
|
|
|
[ 1 z ] [y] = [d]
|
|
[ 0 B ] [x] [a]
|
|
|
|
Using the Inv(B~) expression earlier, the FTRAN result is therefore:
|
|
|
|
[y] = [ 1 -z*Inv(B) ] [d] = [ d - z*Inv(B)*a ]
|
|
[x] [ 0 Inv(B) ] [a] [ Inv(B)*a ]
|
|
|
|
As an example, the value of the dual objective can be returned at the
|
|
0-th index by passing the active RHS vector with 0 at the 0-th position.
|
|
|
|
Similarily, doing the left solve - performing the BTRAN calculation:
|
|
|
|
[x y] [ 1 z ] = [d a']
|
|
[ 0 B ]
|
|
|
|
... will produce the following result in lp_solve:
|
|
|
|
[x y] = [d a'] [ 1 -z*Inv(B) ] = [ d | -d*z*Inv(B) + a'*Inv(B) ]
|
|
[ 0 Inv(B) ]
|
|
|
|
So, if you thought you were simply computing "a'*Inv(B)", look again.
|
|
In order to produce the desired result, you have to set d to 0 before
|
|
the BTRAN operation. On the other hand, if you set d to 1 and a to 0,
|
|
then you are very conveniently on your way to obtain the reduced costs
|
|
(needs a further matrix premultiplication with non-basic variables).
|
|
|
|
Incidentally, the BTRAN with [1 0] that yields [ 1 | -z*Inv(B) ] can
|
|
also be used as a fast way of checking the accuracy of the current
|
|
factorization.
|
|
|
|
Equipped with this understanding, I hope that you see that
|
|
the approach in lp_solve is actually pretty convenient. It also
|
|
becomes easier to extend functionality in lp_solve by drawing on
|
|
formulas and expressions from LP literature that otherwise assume
|
|
the non-bordered syntax and representation.
|
|
|
|
Kjell Eikland -- November 2003
|
|
KE update -- April 2005
|
|
KE update -- June 2005
|
|
|
|
*/
|
|
|
|
STATIC MYBOOL __WINAPI invert(lprec *lp, MYBOOL shiftbounds, MYBOOL final)
|
|
{
|
|
MYBOOL *usedpos, resetbasis;
|
|
REAL test;
|
|
int k, i, j;
|
|
int singularities, usercolB;
|
|
|
|
/* Make sure the tags are correct */
|
|
if(!mat_validate(lp->matA)) {
|
|
lp->spx_status = INFEASIBLE;
|
|
return(FALSE);
|
|
}
|
|
|
|
/* Create the inverse management object at the first call to invert() */
|
|
if(lp->invB == NULL)
|
|
lp->bfp_init(lp, lp->rows, 0, NULL);
|
|
else
|
|
lp->bfp_preparefactorization(lp);
|
|
singularities = 0;
|
|
|
|
/* Must save spx_status since it is used to carry information about
|
|
the presence and handling of singular columns in the matrix */
|
|
if(userabort(lp, MSG_INVERT))
|
|
return(FALSE);
|
|
|
|
#ifdef Paranoia
|
|
if(lp->spx_trace)
|
|
report(lp, DETAILED, "invert: Iter %10g, fact-length %7d, OF " RESULTVALUEMASK ".\n",
|
|
(double) get_total_iter(lp), lp->bfp_colcount(lp), (double) -lp->rhs[0]);
|
|
#endif
|
|
|
|
/* Store state of pre-existing basis, and at the same time check if
|
|
the basis is I; in this case take the easy way out */
|
|
if(!allocMYBOOL(lp, &usedpos, lp->sum + 1, TRUE)) {
|
|
lp->bb_break = TRUE;
|
|
return(FALSE);
|
|
}
|
|
usedpos[0] = TRUE;
|
|
usercolB = 0;
|
|
for(i = 1; i <= lp->rows; i++) {
|
|
k = lp->var_basic[i];
|
|
if(k > lp->rows)
|
|
usercolB++;
|
|
usedpos[k] = TRUE;
|
|
}
|
|
#ifdef Paranoia
|
|
if(!verify_basis(lp))
|
|
report(lp, SEVERE, "invert: Invalid basis detected (iter %g).\n",
|
|
(double) get_total_iter(lp));
|
|
#endif
|
|
|
|
/* Tally matrix nz-counts and check if we should reset basis
|
|
indicators to all slacks */
|
|
resetbasis = (MYBOOL) ((usercolB > 0) && lp->bfp_canresetbasis(lp));
|
|
k = 0;
|
|
for(i = 1; i <= lp->rows; i++) {
|
|
if(lp->var_basic[i] > lp->rows)
|
|
k += mat_collength(lp->matA, lp->var_basic[i] - lp->rows) + (is_OF_nz(lp,lp->var_basic[i] - lp->rows) ? 1 : 0);
|
|
if(resetbasis) {
|
|
j = lp->var_basic[i];
|
|
if(j > lp->rows)
|
|
lp->is_basic[j] = FALSE;
|
|
lp->var_basic[i] = i;
|
|
lp->is_basic[i] = TRUE;
|
|
}
|
|
}
|
|
|
|
/* Now do the refactorization */
|
|
singularities = lp->bfp_factorize(lp, usercolB, k, usedpos, final);
|
|
|
|
/* Do user reporting */
|
|
if(userabort(lp, MSG_INVERT))
|
|
goto Cleanup;
|
|
|
|
/* Finalize factorization/inversion */
|
|
lp->bfp_finishfactorization(lp);
|
|
|
|
/* Recompute the RHS ( Ref. lp_solve inverse logic and Chvatal p. 121 ) */
|
|
#ifdef DebugInv
|
|
blockWriteLREAL(stdout, "RHS-values pre invert", lp->rhs, 0, lp->rows);
|
|
#endif
|
|
recompute_solution(lp, shiftbounds);
|
|
restartPricer(lp, AUTOMATIC);
|
|
#ifdef DebugInv
|
|
blockWriteLREAL(stdout, "RHS-values post invert", lp->rhs, 0, lp->rows);
|
|
#endif
|
|
|
|
Cleanup:
|
|
/* Check for numerical instability indicated by frequent refactorizations */
|
|
test = get_refactfrequency(lp, FALSE);
|
|
if(test < MIN_REFACTFREQUENCY) {
|
|
test = get_refactfrequency(lp, TRUE);
|
|
report(lp, NORMAL, "invert: Refactorization frequency %.1g indicates numeric instability.\n",
|
|
test);
|
|
lp->spx_status = NUMFAILURE;
|
|
}
|
|
|
|
FREE(usedpos);
|
|
return((MYBOOL) (singularities <= 0));
|
|
} /* invert */
|
|
|
|
|
|
STATIC MYBOOL fimprove(lprec *lp, REAL *pcol, int *nzidx, REAL roundzero)
|
|
{
|
|
REAL *errors, sdp;
|
|
int j;
|
|
MYBOOL Ok = TRUE;
|
|
|
|
allocREAL(lp, &errors, lp->rows + 1, FALSE);
|
|
if(errors == NULL) {
|
|
Ok = FALSE;
|
|
return(Ok);
|
|
}
|
|
MEMCOPY(errors, pcol, lp->rows + 1);
|
|
lp->bfp_ftran_normal(lp, pcol, nzidx);
|
|
prod_Ax(lp, NULL, pcol, NULL, 0.0, -1,
|
|
errors, NULL, MAT_ROUNDDEFAULT);
|
|
lp->bfp_ftran_normal(lp, errors, NULL);
|
|
|
|
sdp = 0;
|
|
for(j = 1; j <= lp->rows; j++)
|
|
if(fabs(errors[j])>sdp)
|
|
sdp = fabs(errors[j]);
|
|
if(sdp > lp->epsmachine) {
|
|
report(lp, DETAILED, "Iterative FTRAN correction metric %g", sdp);
|
|
for(j = 1; j <= lp->rows; j++) {
|
|
pcol[j] += errors[j];
|
|
my_roundzero(pcol[j], roundzero);
|
|
}
|
|
}
|
|
FREE(errors);
|
|
return(Ok);
|
|
}
|
|
|
|
STATIC MYBOOL bimprove(lprec *lp, REAL *rhsvector, int *nzidx, REAL roundzero)
|
|
{
|
|
int j;
|
|
REAL *errors, err, maxerr;
|
|
MYBOOL Ok = TRUE;
|
|
|
|
allocREAL(lp, &errors, lp->sum + 1, FALSE);
|
|
if(errors == NULL) {
|
|
Ok = FALSE;
|
|
return(Ok);
|
|
}
|
|
MEMCOPY(errors, rhsvector, lp->sum + 1);
|
|
|
|
/* Solve Ax=b for x, compute b back */
|
|
lp->bfp_btran_normal(lp, errors, nzidx);
|
|
prod_xA(lp, NULL, errors, NULL, 0.0, 1.0,
|
|
errors, NULL,
|
|
MAT_ROUNDDEFAULT);
|
|
|
|
/* Take difference with ingoing values, while shifting the column values
|
|
to the rows section and zeroing the columns again */
|
|
for(j = 1; j <= lp->rows; j++)
|
|
errors[j] = errors[lp->rows+lp->var_basic[j]] - rhsvector[j];
|
|
for(j = lp->rows; j <= lp->sum; j++)
|
|
errors[j] = 0;
|
|
|
|
/* Solve the b errors for the iterative x adjustment */
|
|
lp->bfp_btran_normal(lp, errors, NULL);
|
|
|
|
/* Generate the adjustments and compute statistic */
|
|
maxerr = 0;
|
|
for(j = 1; j <= lp->rows; j++) {
|
|
if(lp->var_basic[j]<=lp->rows) continue;
|
|
err = errors[lp->rows+lp->var_basic[j]];
|
|
if(fabs(err)>maxerr)
|
|
maxerr = fabs(err);
|
|
}
|
|
if(maxerr > lp->epsmachine) {
|
|
report(lp, DETAILED, "Iterative BTRAN correction metric %g", maxerr);
|
|
for(j = 1; j <= lp->rows; j++) {
|
|
if(lp->var_basic[j]<=lp->rows) continue;
|
|
rhsvector[j] += errors[lp->rows+lp->var_basic[j]];
|
|
my_roundzero(rhsvector[j], roundzero);
|
|
}
|
|
}
|
|
FREE(errors);
|
|
return(Ok);
|
|
}
|
|
|
|
STATIC void ftran(lprec *lp, REAL *rhsvector, int *nzidx, REAL roundzero)
|
|
{
|
|
#if 0
|
|
if(is_action(lp->improve, IMPROVE_SOLUTION) && lp->bfp_pivotcount(lp))
|
|
fimprove(lp, rhsvector, nzidx, roundzero);
|
|
else
|
|
#endif
|
|
lp->bfp_ftran_normal(lp, rhsvector, nzidx);
|
|
}
|
|
|
|
STATIC void btran(lprec *lp, REAL *rhsvector, int *nzidx, REAL roundzero)
|
|
{
|
|
#if 0
|
|
if(is_action(lp->improve, IMPROVE_SOLUTION) && lp->bfp_pivotcount(lp))
|
|
bimprove(lp, rhsvector, nzidx, roundzero);
|
|
else
|
|
#endif
|
|
lp->bfp_btran_normal(lp, rhsvector, nzidx);
|
|
}
|
|
|
|
STATIC MYBOOL fsolve(lprec *lp, int varin, REAL *pcol, int *nzidx, REAL roundzero, REAL ofscalar, MYBOOL prepareupdate)
|
|
/* Was setpivcol in versions earlier than 4.0.1.8 - KE */
|
|
{
|
|
MYBOOL ok = TRUE;
|
|
|
|
if(varin > 0)
|
|
obtain_column(lp, varin, pcol, nzidx, NULL);
|
|
|
|
/* Solve, adjusted for objective function scalar */
|
|
pcol[0] *= ofscalar;
|
|
if(prepareupdate)
|
|
lp->bfp_ftran_prepare(lp, pcol, nzidx);
|
|
else
|
|
ftran(lp, pcol, nzidx, roundzero);
|
|
|
|
return(ok);
|
|
|
|
} /* fsolve */
|
|
|
|
|
|
STATIC MYBOOL bsolve(lprec *lp, int row_nr, REAL *rhsvector, int *nzidx, REAL roundzero, REAL ofscalar)
|
|
{
|
|
MYBOOL ok = TRUE;
|
|
|
|
if(row_nr >= 0) /* Note that row_nr == 0 returns the [1, 0...0 ] vector */
|
|
row_nr = obtain_column(lp, row_nr, rhsvector, nzidx, NULL);
|
|
|
|
/* Solve, adjusted for objective function scalar */
|
|
rhsvector[0] *= ofscalar;
|
|
btran(lp, rhsvector, nzidx, roundzero);
|
|
|
|
return(ok);
|
|
|
|
} /* bsolve */
|
|
|
|
|
|
/* Vector compression and expansion routines */
|
|
STATIC MYBOOL vec_compress(REAL *densevector, int startpos, int endpos, REAL epsilon,
|
|
REAL *nzvector, int *nzindex)
|
|
{
|
|
int n;
|
|
|
|
if((densevector == NULL) || (nzindex == NULL) || (startpos > endpos))
|
|
return( FALSE );
|
|
|
|
n = 0;
|
|
densevector += startpos;
|
|
while(startpos <= endpos) {
|
|
if(fabs(*densevector) > epsilon) { /* Apply zero-threshold */
|
|
if(nzvector != NULL) /* Only produce index if no nzvector is given */
|
|
nzvector[n] = *densevector;
|
|
n++;
|
|
nzindex[n] = startpos;
|
|
}
|
|
startpos++;
|
|
densevector++;
|
|
}
|
|
nzindex[0] = n;
|
|
return( TRUE );
|
|
}
|
|
|
|
STATIC MYBOOL vec_expand(REAL *nzvector, int *nzindex, REAL *densevector, int startpos, int endpos)
|
|
{
|
|
int i, n;
|
|
|
|
n = nzindex[0];
|
|
i = nzindex[n];
|
|
densevector += endpos;
|
|
while(endpos >= startpos) { /* Loop from behind to allow densevector == nzvector */
|
|
if(endpos == i) {
|
|
n--;
|
|
*densevector = nzvector[n];
|
|
i = nzindex[n];
|
|
}
|
|
else
|
|
*densevector = 0;
|
|
endpos--;
|
|
densevector--;
|
|
}
|
|
return( TRUE );
|
|
}
|
|
|
|
|
|
/* ----------------------------------------------------------------------- */
|
|
/* Sparse matrix product routines and utility */
|
|
/* ----------------------------------------------------------------------- */
|
|
|
|
STATIC MYBOOL get_colIndexA(lprec *lp, int varset, int *colindex, MYBOOL append)
|
|
{
|
|
int i, varnr, P1extraDim, vb, ve, n, nrows = lp->rows, nsum = lp->sum;
|
|
MYBOOL omitfixed, omitnonfixed;
|
|
REAL v;
|
|
|
|
/* Find what variable range to scan - default is {SCAN_USERVARS} */
|
|
/* First determine the starting position; add from the top, going down */
|
|
P1extraDim = abs(lp->P1extraDim);
|
|
vb = nrows + 1;
|
|
if(varset & SCAN_ARTIFICIALVARS)
|
|
vb = nsum - P1extraDim + 1;
|
|
if(varset & SCAN_USERVARS)
|
|
vb = nrows + 1;
|
|
if(varset & SCAN_SLACKVARS)
|
|
vb = 1;
|
|
|
|
/* Then determine the ending position, add from the bottom, going up */
|
|
ve = nsum;
|
|
if(varset & SCAN_SLACKVARS)
|
|
ve = nrows;
|
|
if(varset & SCAN_USERVARS)
|
|
ve = nsum - P1extraDim;
|
|
if(varset & SCAN_ARTIFICIALVARS)
|
|
ve = nsum;
|
|
|
|
/* Adjust for partial pricing */
|
|
if(varset & SCAN_PARTIALBLOCK) {
|
|
SETMAX(vb, partial_blockStart(lp, FALSE));
|
|
SETMIN(ve, partial_blockEnd(lp, FALSE));
|
|
}
|
|
|
|
/* Determine exclusion columns */
|
|
omitfixed = (MYBOOL) ((varset & OMIT_FIXED) != 0);
|
|
omitnonfixed = (MYBOOL) ((varset & OMIT_NONFIXED) != 0);
|
|
if(omitfixed && omitnonfixed)
|
|
return(FALSE);
|
|
|
|
/* Scan the target colums */
|
|
if(append)
|
|
n = colindex[0];
|
|
else
|
|
n = 0;
|
|
for(varnr = vb; varnr <= ve; varnr++) {
|
|
|
|
/* Skip gap in the specified column scan range (possibly user variables) */
|
|
if(varnr > nrows) {
|
|
if((varnr <= nsum-P1extraDim) && !(varset & SCAN_USERVARS))
|
|
continue;
|
|
#if 1
|
|
/* Skip empty columns */
|
|
if(/*(lp->P1extraVal == 0) &&*/
|
|
(mat_collength(lp->matA, varnr-nrows) == 0))
|
|
continue;
|
|
#endif
|
|
}
|
|
|
|
/* Find if the variable is in the scope - default is {�} */
|
|
i = lp->is_basic[varnr];
|
|
if((varset & USE_BASICVARS) > 0 && (i))
|
|
;
|
|
else if((varset & USE_NONBASICVARS) > 0 && (!i))
|
|
;
|
|
else
|
|
continue;
|
|
|
|
v = lp->upbo[varnr];
|
|
if((omitfixed && (v == 0)) ||
|
|
(omitnonfixed && (v != 0)))
|
|
continue;
|
|
|
|
/* Append to list */
|
|
n++;
|
|
colindex[n] = varnr;
|
|
}
|
|
colindex[0] = n;
|
|
|
|
return(TRUE);
|
|
}
|
|
|
|
STATIC int prod_Ax(lprec *lp, int *coltarget, REAL *input, int *nzinput,
|
|
REAL roundzero, REAL ofscalar,
|
|
REAL *output, int *nzoutput, int roundmode)
|
|
/* prod_Ax is only used in fimprove; note that it is NOT VALIDATED/verified as of 20030801 - KE */
|
|
{
|
|
int j, colnr, ib, ie, vb, ve;
|
|
MYBOOL localset, localnz = FALSE, isRC;
|
|
MATrec *mat = lp->matA;
|
|
REAL sdp;
|
|
REAL *value;
|
|
int *rownr;
|
|
|
|
/* Find what variable range to scan - default is {SCAN_USERVARS} */
|
|
/* Define default column target if none was provided */
|
|
isRC = (MYBOOL) ((roundmode & MAT_ROUNDRC) != 0);
|
|
localset = (MYBOOL) (coltarget == NULL);
|
|
if(localset) {
|
|
int varset = SCAN_SLACKVARS | SCAN_USERVARS |
|
|
USE_BASICVARS | OMIT_FIXED;
|
|
if(isRC && is_piv_mode(lp, PRICE_PARTIAL) && !is_piv_mode(lp, PRICE_FORCEFULL))
|
|
varset |= SCAN_PARTIALBLOCK;
|
|
coltarget = (int *) mempool_obtainVector(lp->workarrays, lp->sum+1, sizeof(*coltarget));
|
|
if(!get_colIndexA(lp, varset, coltarget, FALSE)) {
|
|
mempool_releaseVector(lp->workarrays, (char *) coltarget, FALSE);
|
|
return(FALSE);
|
|
}
|
|
}
|
|
localnz = (MYBOOL) (nzinput == NULL);
|
|
if(localnz) {
|
|
nzinput = (int *) mempool_obtainVector(lp->workarrays, lp->rows+1, sizeof(*nzinput));
|
|
vec_compress(input, 0, lp->rows, lp->matA->epsvalue, NULL, nzinput);
|
|
}
|
|
|
|
/* Scan the columns */
|
|
vb = 1;
|
|
ve = coltarget[0];
|
|
for(vb = 1; vb <= coltarget[0]; vb++) {
|
|
colnr = coltarget[vb];
|
|
j = lp->is_basic[colnr];
|
|
|
|
/* Perform the multiplication */
|
|
sdp = ofscalar*input[j];
|
|
if(colnr <= lp->rows) /* A slack variable is in the basis */
|
|
output[colnr] += sdp;
|
|
else { /* A normal variable is in the basis */
|
|
colnr -= lp->rows;
|
|
ib = mat->col_end[colnr - 1];
|
|
ie = mat->col_end[colnr];
|
|
rownr = &COL_MAT_ROWNR(ib);
|
|
value = &COL_MAT_VALUE(ib);
|
|
for(; ib < ie;
|
|
ib++, rownr += matRowColStep, value += matValueStep) {
|
|
output[*rownr] += (*value)*sdp;
|
|
}
|
|
}
|
|
}
|
|
roundVector(output+1, lp->rows-1, roundzero);
|
|
|
|
/* Clean up and return */
|
|
if(localset)
|
|
mempool_releaseVector(lp->workarrays, (char *) coltarget, FALSE);
|
|
if(localnz)
|
|
mempool_releaseVector(lp->workarrays, (char *) nzinput, FALSE);
|
|
|
|
return(TRUE);
|
|
}
|
|
|
|
STATIC int prod_xA(lprec *lp, int *coltarget,
|
|
REAL *input, int *nzinput, REAL roundzero, REAL ofscalar,
|
|
REAL *output, int *nzoutput, int roundmode)
|
|
/* Note that the dot product xa is stored at the active column index of A, i.e. of a.
|
|
This means that if the basis only contains non-slack variables, output may point to
|
|
the same vector as input, without overwriting the [0..rows] elements. */
|
|
{
|
|
int colnr, rownr, varnr, ib, ie, vb, ve, nrows = lp->rows;
|
|
MYBOOL localset, localnz = FALSE, includeOF, isRC;
|
|
REALXP vmax;
|
|
register REALXP v;
|
|
int inz, *rowin, countNZ = 0;
|
|
MATrec *mat = lp->matA;
|
|
register REAL *matValue;
|
|
register int *matRownr;
|
|
|
|
/* Clean output area (only necessary if we are returning the full vector) */
|
|
isRC = (MYBOOL) ((roundmode & MAT_ROUNDRC) != 0);
|
|
if(nzoutput == NULL) {
|
|
if(input == output)
|
|
MEMCLEAR(output+nrows+1, lp->columns);
|
|
else
|
|
MEMCLEAR(output, lp->sum+1);
|
|
}
|
|
|
|
/* Find what variable range to scan - default is {SCAN_USERVARS} */
|
|
/* Define default column target if none was provided */
|
|
localset = (MYBOOL) (coltarget == NULL);
|
|
if(localset) {
|
|
int varset = SCAN_SLACKVARS | SCAN_USERVARS |
|
|
USE_NONBASICVARS | OMIT_FIXED;
|
|
if(isRC && is_piv_mode(lp, PRICE_PARTIAL) && !is_piv_mode(lp, PRICE_FORCEFULL))
|
|
varset |= SCAN_PARTIALBLOCK;
|
|
coltarget = (int *) mempool_obtainVector(lp->workarrays, lp->sum+1, sizeof(*coltarget));
|
|
if(!get_colIndexA(lp, varset, coltarget, FALSE)) {
|
|
mempool_releaseVector(lp->workarrays, (char *) coltarget, FALSE);
|
|
return(FALSE);
|
|
}
|
|
}
|
|
/*#define UseLocalNZ*/
|
|
#ifdef UseLocalNZ
|
|
localnz = (MYBOOL) (nzinput == NULL);
|
|
if(localnz) {
|
|
nzinput = (int *) mempool_obtainVector(lp->workarrays, nrows+1, sizeof(*nzinput));
|
|
vec_compress(input, 0, nrows, lp->matA->epsvalue, NULL, nzinput);
|
|
}
|
|
#endif
|
|
includeOF = (MYBOOL) (((nzinput == NULL) || (nzinput[1] == 0)) &&
|
|
(input[0] != 0) && lp->obj_in_basis);
|
|
|
|
/* Scan the target colums */
|
|
vmax = 0;
|
|
ve = coltarget[0];
|
|
for(vb = 1; vb <= ve; vb++) {
|
|
|
|
varnr = coltarget[vb];
|
|
|
|
if(varnr <= nrows) {
|
|
v = input[varnr];
|
|
}
|
|
else {
|
|
colnr = varnr - nrows;
|
|
v = 0;
|
|
ib = mat->col_end[colnr - 1];
|
|
ie = mat->col_end[colnr];
|
|
if(ib < ie) {
|
|
|
|
/* Do dense input vector version */
|
|
#ifdef UseLocalNZ
|
|
if(localnz || (nzinput == NULL)) {
|
|
#else
|
|
if(nzinput == NULL) {
|
|
#endif
|
|
/* Do the OF */
|
|
if(includeOF)
|
|
#ifdef DirectArrayOF
|
|
v += input[0] * lp->obj[colnr] * ofscalar;
|
|
#else
|
|
v += input[0] * get_OF_active(lp, varnr, ofscalar);
|
|
#endif
|
|
|
|
/* Initialize pointers */
|
|
matRownr = &COL_MAT_ROWNR(ib);
|
|
matValue = &COL_MAT_VALUE(ib);
|
|
|
|
/* Do extra loop optimization based on target window overlaps */
|
|
#ifdef UseLocalNZ
|
|
if((ib < ie)
|
|
&& (colnr <= *nzinput)
|
|
&& (COL_MAT_ROWNR(ie-1) >= nzinput[colnr])
|
|
&& (*matRownr <= nzinput[*nzinput])
|
|
)
|
|
#endif
|
|
#ifdef NoLoopUnroll
|
|
/* Then loop over all regular rows */
|
|
for(; ib < ie; ib++) {
|
|
v += input[*matRownr] * (*matValue);
|
|
matValue += matValueStep;
|
|
matRownr += matRowColStep;
|
|
}
|
|
#else
|
|
/* Prepare for simple loop unrolling */
|
|
if(((ie-ib) % 2) == 1) {
|
|
v += input[*matRownr] * (*matValue);
|
|
ib++;
|
|
matValue += matValueStep;
|
|
matRownr += matRowColStep;
|
|
}
|
|
|
|
/* Then loop over remaining pairs of regular rows */
|
|
while(ib < ie) {
|
|
v += input[*matRownr] * (*matValue);
|
|
v += input[*(matRownr+matRowColStep)] * (*(matValue+matValueStep));
|
|
ib += 2;
|
|
matValue += 2*matValueStep;
|
|
matRownr += 2*matRowColStep;
|
|
}
|
|
#endif
|
|
}
|
|
/* Do sparse input vector version */
|
|
else {
|
|
|
|
/* Do the OF */
|
|
if(includeOF)
|
|
#ifdef DirectArrayOF
|
|
v += input[0] * lp->obj[colnr] * ofscalar;
|
|
#else
|
|
v += input[0] * get_OF_active(lp, varnr, ofscalar);
|
|
#endif
|
|
|
|
/* Initialize pointers */
|
|
inz = 1;
|
|
rowin = nzinput+inz;
|
|
matRownr = &COL_MAT_ROWNR(ib);
|
|
matValue = &COL_MAT_VALUE(ib);
|
|
ie--;
|
|
|
|
/* Then loop over all non-OF rows */
|
|
while((inz <= *nzinput) && (ib <= ie)) {
|
|
|
|
/* Try to synchronize at right */
|
|
while((*rowin > *matRownr) && (ib < ie)) {
|
|
ib++;
|
|
matValue += matValueStep;
|
|
matRownr += matRowColStep;
|
|
}
|
|
/* Try to synchronize at left */
|
|
while((*rowin < *matRownr) && (inz < *nzinput)) {
|
|
inz++;
|
|
rowin++;
|
|
}
|
|
/* Perform dot product operation if there was a match */
|
|
if(*rowin == *matRownr) {
|
|
v += input[*rowin] * (*matValue);
|
|
/* Step forward at left */
|
|
inz++;
|
|
rowin++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
if((roundmode & MAT_ROUNDABS) != 0) {
|
|
my_roundzero(v, roundzero);
|
|
}
|
|
}
|
|
|
|
/* Special handling of small reduced cost values */
|
|
if(!isRC || (my_chsign(lp->is_lower[varnr], v) < 0)) {
|
|
SETMAX(vmax, fabs((REAL) v));
|
|
}
|
|
if(v != 0) {
|
|
countNZ++;
|
|
if(nzoutput != NULL)
|
|
nzoutput[countNZ] = varnr;
|
|
}
|
|
output[varnr] = (REAL) v;
|
|
}
|
|
|
|
/* Compute reduced cost if this option is active */
|
|
if(isRC && !lp->obj_in_basis)
|
|
countNZ = get_basisOF(lp, coltarget, output, nzoutput);
|
|
|
|
/* Check if we should do relative rounding */
|
|
if((roundmode & MAT_ROUNDREL) != 0) {
|
|
if((roundzero > 0) && (nzoutput != NULL)) {
|
|
ie = 0;
|
|
if(isRC) {
|
|
SETMAX(vmax, MAT_ROUNDRCMIN); /* Make sure we don't use very small values */
|
|
}
|
|
vmax *= roundzero;
|
|
for(ib = 1; ib <= countNZ; ib++) {
|
|
rownr = nzoutput[ib];
|
|
if(fabs(output[rownr]) < vmax)
|
|
output[rownr] = 0;
|
|
else {
|
|
ie++;
|
|
nzoutput[ie] = rownr;
|
|
}
|
|
}
|
|
countNZ = ie;
|
|
}
|
|
}
|
|
|
|
/* Clean up and return */
|
|
if(localset)
|
|
mempool_releaseVector(lp->workarrays, (char *) coltarget, FALSE);
|
|
if(localnz)
|
|
mempool_releaseVector(lp->workarrays, (char *) nzinput, FALSE);
|
|
|
|
if(nzoutput != NULL)
|
|
*nzoutput = countNZ;
|
|
return(countNZ);
|
|
}
|
|
|
|
STATIC MYBOOL prod_xA2(lprec *lp, int *coltarget,
|
|
REAL *prow, REAL proundzero, int *nzprow,
|
|
REAL *drow, REAL droundzero, int *nzdrow,
|
|
REAL ofscalar, int roundmode)
|
|
{
|
|
int varnr, colnr, ib, ie, vb, ve, nrows = lp->rows;
|
|
MYBOOL includeOF, isRC;
|
|
REALXP dmax, pmax;
|
|
register REALXP d, p;
|
|
MATrec *mat = lp->matA;
|
|
REAL value;
|
|
register REAL *matValue;
|
|
register int *matRownr;
|
|
MYBOOL localset;
|
|
|
|
/* Find what variable range to scan - default is {SCAN_USERVARS} */
|
|
/* First determine the starting position; add from the top, going down */
|
|
localset = (MYBOOL) (coltarget == NULL);
|
|
if(localset) {
|
|
int varset = SCAN_SLACKVARS + SCAN_USERVARS + /*SCAN_ALLVARS +*/
|
|
/*SCAN_PARTIALBLOCK+*/
|
|
USE_NONBASICVARS+OMIT_FIXED;
|
|
coltarget = (int *) mempool_obtainVector(lp->workarrays, lp->sum+1, sizeof(*coltarget));
|
|
if(!get_colIndexA(lp, varset, coltarget, FALSE)) {
|
|
mempool_releaseVector(lp->workarrays, (char *) coltarget, FALSE);
|
|
return(FALSE);
|
|
}
|
|
}
|
|
|
|
/* Initialize variables */
|
|
isRC = (MYBOOL) ((roundmode & MAT_ROUNDRC) != 0);
|
|
pmax = 0;
|
|
dmax = 0;
|
|
if(nzprow != NULL)
|
|
*nzprow = 0;
|
|
if(nzdrow != NULL)
|
|
*nzdrow = 0;
|
|
includeOF = (MYBOOL) (((prow[0] != 0) || (drow[0] != 0)) &&
|
|
lp->obj_in_basis);
|
|
|
|
/* Scan the target colums */
|
|
ve = coltarget[0];
|
|
for(vb = 1; vb <= ve; vb++) {
|
|
|
|
varnr = coltarget[vb];
|
|
|
|
if(varnr <= nrows) {
|
|
p = prow[varnr];
|
|
d = drow[varnr];
|
|
}
|
|
else {
|
|
|
|
colnr = varnr - nrows;
|
|
|
|
p = 0;
|
|
d = 0;
|
|
ib = mat->col_end[colnr - 1];
|
|
ie = mat->col_end[colnr];
|
|
|
|
if(ib < ie) {
|
|
|
|
/* Do the OF */
|
|
if(includeOF) {
|
|
#ifdef DirectArrayOF
|
|
value = lp->obj[colnr] * ofscalar;
|
|
#else
|
|
value = get_OF_active(lp, varnr, ofscalar);
|
|
#endif
|
|
p += prow[0] * value;
|
|
d += drow[0] * value;
|
|
}
|
|
|
|
/* Then loop over all regular rows */
|
|
matRownr = &COL_MAT_ROWNR(ib);
|
|
matValue = &COL_MAT_VALUE(ib);
|
|
#ifdef NoLoopUnroll
|
|
for( ; ib < ie; ib++) {
|
|
p += prow[*matRownr] * (*matValue);
|
|
d += drow[*matRownr] * (*matValue);
|
|
matValue += matValueStep;
|
|
matRownr += matRowColStep;
|
|
}
|
|
#else
|
|
/* Prepare for simple loop unrolling */
|
|
if(((ie-ib) % 2) == 1) {
|
|
p += prow[*matRownr] * (*matValue);
|
|
d += drow[*matRownr] * (*matValue);
|
|
ib++;
|
|
matValue += matValueStep;
|
|
matRownr += matRowColStep;
|
|
}
|
|
|
|
/* Then loop over remaining pairs of regular rows */
|
|
while(ib < ie) {
|
|
p += prow[*matRownr] * (*matValue);
|
|
p += prow[*(matRownr+matRowColStep)] * (*(matValue+matValueStep));
|
|
d += drow[*matRownr] * (*matValue);
|
|
d += drow[*(matRownr+matRowColStep)] * (*(matValue+matValueStep));
|
|
ib += 2;
|
|
matValue += 2*matValueStep;
|
|
matRownr += 2*matRowColStep;
|
|
}
|
|
#endif
|
|
|
|
}
|
|
if((roundmode & MAT_ROUNDABS) != 0) {
|
|
my_roundzero(p, proundzero);
|
|
my_roundzero(d, droundzero);
|
|
}
|
|
}
|
|
|
|
SETMAX(pmax, fabs((REAL) p));
|
|
prow[varnr] = (REAL) p;
|
|
if((nzprow != NULL) && (p != 0)) {
|
|
(*nzprow)++;
|
|
nzprow[*nzprow] = varnr;
|
|
}
|
|
|
|
/* Special handling of reduced cost rounding */
|
|
if(!isRC || (my_chsign(lp->is_lower[varnr], d) < 0)) {
|
|
SETMAX(dmax, fabs((REAL) d));
|
|
}
|
|
drow[varnr] = (REAL) d;
|
|
if((nzdrow != NULL) && (d != 0)) {
|
|
(*nzdrow)++;
|
|
nzdrow[*nzdrow] = varnr;
|
|
}
|
|
}
|
|
|
|
/* Compute reduced cost here if this option is active */
|
|
if((drow != 0) && !lp->obj_in_basis)
|
|
get_basisOF(lp, coltarget, drow, nzdrow);
|
|
|
|
/* Check if we should do relative rounding */
|
|
if((roundmode & MAT_ROUNDREL) != 0) {
|
|
if((proundzero > 0) && (nzprow != NULL)) {
|
|
ie = 0;
|
|
pmax *= proundzero;
|
|
for(ib = 1; ib <= *nzprow; ib++) {
|
|
varnr = nzprow[ib];
|
|
if(fabs(prow[varnr]) < pmax)
|
|
prow[varnr] = 0;
|
|
else {
|
|
ie++;
|
|
nzprow[ie] = varnr;
|
|
}
|
|
}
|
|
*nzprow = ie;
|
|
}
|
|
if((droundzero > 0) && (nzdrow != NULL)) {
|
|
ie = 0;
|
|
if(isRC) {
|
|
SETMAX(dmax, MAT_ROUNDRCMIN); /* Make sure we don't use very small values */
|
|
}
|
|
dmax *= droundzero;
|
|
for(ib = 1; ib <= *nzdrow; ib++) {
|
|
varnr = nzdrow[ib];
|
|
if(fabs(drow[varnr]) < dmax)
|
|
drow[varnr] = 0;
|
|
else {
|
|
ie++;
|
|
nzdrow[ie] = varnr;
|
|
}
|
|
}
|
|
*nzdrow = ie;
|
|
}
|
|
}
|
|
|
|
/* Clean up and return */
|
|
if(localset)
|
|
mempool_releaseVector(lp->workarrays, (char *) coltarget, FALSE);
|
|
return( TRUE );
|
|
}
|
|
|
|
STATIC void bsolve_xA2(lprec *lp, int* coltarget,
|
|
int row_nr1, REAL *vector1, REAL roundzero1, int *nzvector1,
|
|
int row_nr2, REAL *vector2, REAL roundzero2, int *nzvector2, int roundmode)
|
|
{
|
|
REAL ofscalar = 1.0;
|
|
|
|
/* Clear and initialize first vector */
|
|
if(nzvector1 == NULL)
|
|
MEMCLEAR(vector1, lp->sum + 1);
|
|
else
|
|
MEMCLEAR(vector1, lp->rows + 1);
|
|
vector1[row_nr1] = 1;
|
|
/* workINT[0] = 1;
|
|
workINT[1] = row_nr1; */
|
|
|
|
if(vector2 == NULL) {
|
|
lp->bfp_btran_normal(lp, vector1, NULL);
|
|
prod_xA(lp, coltarget, vector1, NULL, roundzero1, ofscalar*0,
|
|
vector1, nzvector1, roundmode);
|
|
}
|
|
else {
|
|
|
|
/* Clear and initialize second vector */
|
|
if(nzvector2 == NULL)
|
|
MEMCLEAR(vector2, lp->sum + 1);
|
|
else
|
|
MEMCLEAR(vector2, lp->rows + 1);
|
|
if(lp->obj_in_basis || (row_nr2 > 0)) {
|
|
vector2[row_nr2] = 1;
|
|
/* workINT[2] = 1;
|
|
workINT[3] = row_nr2; */
|
|
}
|
|
else
|
|
get_basisOF(lp, NULL, vector2, nzvector2);
|
|
|
|
/* A double BTRAN equation solver process is implemented "in-line" below in
|
|
order to save time and to implement different rounding for the two */
|
|
lp->bfp_btran_double(lp, vector1, NULL, vector2, NULL);
|
|
|
|
/* Multiply solution vectors with matrix values */
|
|
prod_xA2(lp, coltarget, vector1, roundzero1, nzvector1,
|
|
vector2, roundzero2, nzvector2,
|
|
ofscalar, roundmode);
|
|
}
|
|
}
|
|
|