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.
241 lines
6.8 KiB
241 lines
6.8 KiB
/*
|
|
Minimum matrix inverse fill-in modules - interface for lp_solve v5.0+
|
|
----------------------------------------------------------------------------------
|
|
Author: Kjell Eikland
|
|
Contact: kjell.eikland@broadpark.no
|
|
License terms: LGPL.
|
|
|
|
Requires: string.h, colamd.h, lp_lib.h
|
|
|
|
Release notes:
|
|
v1.0 1 September 2003 Preprocessing routines for minimum fill-in column
|
|
ordering for inverse factorization using the open
|
|
source COLAMD library. Suitable for the dense parts
|
|
of both the product form and LU factorization inverse
|
|
methods.
|
|
v1.1 1 July 2004 Renamed from lp_colamdMDO to lp_MDO.
|
|
|
|
----------------------------------------------------------------------------------
|
|
*/
|
|
|
|
#include <string.h>
|
|
#include "commonlib.h"
|
|
#include "lp_lib.h"
|
|
#include "colamd.h"
|
|
#include "lp_MDO.h"
|
|
|
|
#ifdef FORTIFY
|
|
# include "lp_fortify.h"
|
|
#endif
|
|
|
|
STATIC MYBOOL includeMDO(MYBOOL *usedpos, int item)
|
|
{
|
|
/* Legend: TRUE => A basic slack variable already in the basis
|
|
FALSE => A column free for being pivoted in
|
|
AUTOMATIC+TRUE => A row-singleton user column pivoted into the basis
|
|
AUTOMATIC+FALSE => A column-singleton user column pivoted into the basis */
|
|
|
|
/* Handle case where we are processing all columns */
|
|
if(usedpos == NULL)
|
|
return( TRUE );
|
|
|
|
else {
|
|
/* Otherwise do the selective case */
|
|
MYBOOL test = usedpos[item];
|
|
#if 1
|
|
return( test != TRUE );
|
|
#else
|
|
test = test & TRUE;
|
|
return( test == FALSE );
|
|
#endif
|
|
}
|
|
}
|
|
|
|
STATIC int prepareMDO(lprec *lp, MYBOOL *usedpos, int *colorder, int *data, int *rowmap)
|
|
/* This routine prepares data structures for colamd(). It is called twice, the first
|
|
time to count applicable non-zero elements by column, and the second time to fill in
|
|
the row indexes of the non-zero values from the first call. Note that the colamd()
|
|
row index base is 0 (which suits lp_solve fine). */
|
|
{
|
|
int i, ii, j, k, kk;
|
|
int nrows = lp->rows+1, ncols = colorder[0];
|
|
int offset = 0, Bnz = 0, Tnz;
|
|
MYBOOL dotally = (MYBOOL) (rowmap == NULL);
|
|
MATrec *mat = lp->matA;
|
|
REAL hold;
|
|
REAL *value;
|
|
int *rownr;
|
|
|
|
if(dotally)
|
|
data[0] = 0;
|
|
|
|
Tnz = nrows - ncols;
|
|
for(j = 1; j <= ncols; j++) {
|
|
kk = colorder[j];
|
|
|
|
/* Process slacks */
|
|
if(kk <= lp->rows) {
|
|
if(includeMDO(usedpos, kk)) {
|
|
if(!dotally)
|
|
data[Bnz] = rowmap[kk]+offset;
|
|
Bnz++;
|
|
}
|
|
Tnz++;
|
|
}
|
|
/* Process user columns */
|
|
else {
|
|
k = kk - lp->rows;
|
|
i = mat->col_end[k-1];
|
|
ii= mat->col_end[k];
|
|
Tnz += ii-i;
|
|
#ifdef Paranoia
|
|
if(i >= ii)
|
|
lp->report(lp, SEVERE, "prepareMDO: Encountered empty basic column %d\n", k);
|
|
#endif
|
|
|
|
/* Detect if we need to do phase 1 adjustments of zero-valued OF variable */
|
|
rownr = &COL_MAT_ROWNR(i);
|
|
value = &COL_MAT_VALUE(i);
|
|
hold = 0;
|
|
if((*rownr > 0) && includeMDO(usedpos, 0) && modifyOF1(lp, kk, &hold, 1.0)) {
|
|
if(!dotally)
|
|
data[Bnz] = offset;
|
|
Bnz++;
|
|
}
|
|
/* Loop over all NZ-variables */
|
|
for(; i < ii;
|
|
i++, value += matValueStep, rownr += matRowColStep) {
|
|
if(!includeMDO(usedpos, *rownr))
|
|
continue;
|
|
/* See if we need to change phase 1 OF value */
|
|
if(*rownr == 0) {
|
|
hold = *value;
|
|
if(!modifyOF1(lp, kk, &hold, 1.0))
|
|
continue;
|
|
}
|
|
/* Tally uneliminated constraint row values */
|
|
if(!dotally)
|
|
data[Bnz] = rowmap[*rownr]+offset;
|
|
Bnz++;
|
|
}
|
|
}
|
|
if(dotally)
|
|
data[j] = Bnz;
|
|
}
|
|
return( Tnz );
|
|
}
|
|
|
|
STATIC MYBOOL verifyMDO(lprec *lp, int *col_end, int *row_nr, int rowmax, int colmax)
|
|
{
|
|
int i, j, n, err = 0;
|
|
|
|
for(i = 1; i <= colmax; i++) {
|
|
n = 0;
|
|
for(j = col_end[i-1]; (j < col_end[i]) && (err == 0); j++, n++) {
|
|
if(row_nr[j] < 0 || row_nr[j] > rowmax)
|
|
err = 1;
|
|
if(n > 0 && row_nr[j] <= row_nr[j-1])
|
|
err = 2;
|
|
n++;
|
|
}
|
|
}
|
|
if(err != 0)
|
|
lp->report(lp, SEVERE, "verifyMDO: Invalid MDO input structure generated (error %d)\n", err);
|
|
return( (MYBOOL) (err == 0) );
|
|
}
|
|
|
|
void *mdo_calloc(size_t size, size_t count)
|
|
{
|
|
return ( calloc(size, count) );
|
|
}
|
|
void mdo_free(void *mem)
|
|
{
|
|
free( mem );
|
|
}
|
|
|
|
|
|
int __WINAPI getMDO(lprec *lp, MYBOOL *usedpos, int *colorder, int *size, MYBOOL symmetric)
|
|
{
|
|
int error = FALSE;
|
|
int nrows = lp->rows+1, ncols = colorder[0];
|
|
int i, j, kk, n;
|
|
int *col_end, *row_map = NULL;
|
|
int Bnz, Blen, *Brows = NULL;
|
|
int stats[COLAMD_STATS];
|
|
double knobs[COLAMD_KNOBS];
|
|
|
|
/* Tally the non-zero counts of the unused columns/rows of the
|
|
basis matrix and store corresponding "net" starting positions */
|
|
allocINT(lp, &col_end, ncols+1, FALSE);
|
|
n = prepareMDO(lp, usedpos, colorder, col_end, NULL);
|
|
Bnz = col_end[ncols];
|
|
|
|
/* Check that we have unused basic columns, otherwise skip analysis */
|
|
if(ncols == 0 || Bnz == 0)
|
|
goto Transfer;
|
|
|
|
/* Get net number of rows and fill mapper */
|
|
allocINT(lp, &row_map, nrows, FALSE);
|
|
nrows = 0;
|
|
for(i = 0; i <= lp->rows; i++) {
|
|
row_map[i] = i-nrows;
|
|
/* Increment eliminated row counter if necessary */
|
|
if(!includeMDO(usedpos, i))
|
|
nrows++;
|
|
}
|
|
nrows = lp->rows+1 - nrows;
|
|
|
|
/* Store row indeces of non-zero values in the basic columns */
|
|
Blen = colamd_recommended(Bnz, nrows, ncols);
|
|
allocINT(lp, &Brows, Blen, FALSE);
|
|
prepareMDO(lp, usedpos, colorder, Brows, row_map);
|
|
#ifdef Paranoia
|
|
verifyMDO(lp, col_end, Brows, nrows, ncols);
|
|
#endif
|
|
|
|
/* Compute the MDO */
|
|
#if 1
|
|
colamd_set_defaults(knobs);
|
|
knobs [COLAMD_DENSE_ROW] = 0.2+0.2 ; /* default changed for UMFPACK */
|
|
knobs [COLAMD_DENSE_COL] = knobs [COLAMD_DENSE_ROW];
|
|
if(symmetric && (nrows == ncols)) {
|
|
MEMCOPY(colorder, Brows, ncols + 1);
|
|
error = !symamd(nrows, colorder, col_end, Brows, knobs, stats, mdo_calloc, mdo_free);
|
|
}
|
|
else
|
|
error = !colamd(nrows, ncols, Blen, Brows, col_end, knobs, stats);
|
|
#else
|
|
if(symmetric && (nrows == ncols)) {
|
|
MEMCOPY(colorder, Brows, ncols + 1);
|
|
error = !symamd(nrows, colorder, col_end, Brows, knobs, stats, mdo_calloc, mdo_free);
|
|
}
|
|
else
|
|
error = !colamd(nrows, ncols, Blen, Brows, col_end, (double *) NULL, stats);
|
|
#endif
|
|
|
|
/* Transfer the estimated optimal ordering, adjusting for index offsets */
|
|
Transfer:
|
|
if(error)
|
|
error = stats[COLAMD_STATUS];
|
|
else {
|
|
MEMCOPY(Brows, colorder, ncols + 1);
|
|
for(j = 0; j < ncols; j++) {
|
|
kk = col_end[j];
|
|
n = Brows[kk+1];
|
|
colorder[j+1] = n;
|
|
}
|
|
}
|
|
|
|
/* Free temporary vectors */
|
|
FREE(col_end);
|
|
if(row_map != NULL)
|
|
FREE(row_map);
|
|
if(Brows != NULL)
|
|
FREE(Brows);
|
|
|
|
if(size != NULL)
|
|
*size = ncols;
|
|
return( error );
|
|
}
|
|
|
|
|