479 lines
15 KiB
C
479 lines
15 KiB
C
/*************************************************************************
|
|
|
|
This code is derived from cs2 v3.7
|
|
Written by Andrew V. Goldberg and Boris Cherkassky
|
|
Modifications for use in snaphu by Curtis W. Chen
|
|
|
|
Parser for cs2 minimum cost flow solver. Originally written to read
|
|
DIMACS format (text) input files. Modified to parse passed data
|
|
from snaphu. This file is included with a #include from
|
|
snaphu_cs2.c.
|
|
|
|
The cs2 code is used here with permission for strictly noncommerical
|
|
use. The original cs2 source code can be downloaded from
|
|
|
|
http://www.igsystems.com/cs2
|
|
|
|
The original cs2 copyright is stated as follows:
|
|
|
|
COPYRIGHT C 1995 IG Systems, Inc. Permission to use for
|
|
evaluation purposes is granted provided that proper
|
|
acknowledgments are given. For a commercial licence, contact
|
|
igsys@eclipse.net.
|
|
|
|
This software comes with NO WARRANTY, expressed or implied. By way
|
|
of example, but not limitation, we make no representations of
|
|
warranties of merchantability or fitness for any particular
|
|
purpose or that the use of the software components or
|
|
documentation will not infringe any patents, copyrights,
|
|
trademarks, or other rights.
|
|
|
|
Copyright 2002 Board of Trustees, Leland Stanford Jr. University
|
|
|
|
*************************************************************************/
|
|
|
|
|
|
|
|
int cs2mcfparse(residue, rowcost, colcost, nNrow, nNcol, n_ad, m_ad, nodes_ad,
|
|
arcs_ad, node_min_ad, m_c_ad, cap_ad )
|
|
|
|
/* parameters passed to set up network */
|
|
signed char **residue; /* 2D array of residues */
|
|
short **rowcost; /* 2D array of row arc costs */
|
|
short **colcost; /* 2D array of col arc costs */
|
|
long nNrow; /* number of nodes per row */
|
|
long nNcol; /* number of nodes per column */
|
|
|
|
/* these parameters are output */
|
|
long *n_ad; /* address of the number of nodes */
|
|
long *m_ad; /* address of the number of arcs */
|
|
node **nodes_ad; /* address of the array of nodes */
|
|
arc **arcs_ad; /* address of the array of arcs */
|
|
long *node_min_ad; /* address of the minimal node */
|
|
double *m_c_ad; /* maximal arc cost */
|
|
short **cap_ad; /* array of capacities (changed to short) */
|
|
|
|
{
|
|
|
|
|
|
#define ABS( x ) ( (x) >= 0 ) ? (x) : -(x)
|
|
|
|
/* variables added for unwrapping parse */
|
|
unsigned int row, col, dir;
|
|
unsigned long narcs, nnodes, nodectr, arcctr, nresidues;
|
|
long cumsupply, temp;
|
|
|
|
|
|
long inf_cap = 0;
|
|
long n, /* internal number of nodes */
|
|
node_min, /* minimal no of node */
|
|
node_max, /* maximal no of nodes */
|
|
*arc_first, /* internal array for holding
|
|
- node degree
|
|
- position of the first outgoing arc */
|
|
*arc_tail, /* internal array: tails of the arcs */
|
|
/* temporary variables carrying no of nodes */
|
|
head, tail, i;
|
|
|
|
long m, /* internal number of arcs */
|
|
/* temporary variables carrying no of arcs */
|
|
last, arc_num, arc_new_num;
|
|
|
|
node *nodes, /* pointers to the node structure */
|
|
*head_p,
|
|
*ndp,
|
|
*in,
|
|
*jn;
|
|
|
|
arc *arcs, /* pointers to the arc structure */
|
|
*arc_current,
|
|
*arc_new,
|
|
*arc_tmp;
|
|
|
|
long excess, /* supply/demand of the node */
|
|
low, /* lowest flow through the arc */
|
|
acap; /* capacity */
|
|
|
|
long cost; /* arc cost */
|
|
|
|
|
|
double dcost, /* arc cost in double mode */
|
|
m_c; /* maximal arc cost */
|
|
|
|
short *cap; /* array of capacities (changed to short) */
|
|
|
|
double total_p, /* total supply */
|
|
total_n, /* total demand */
|
|
cap_out, /* sum of outgoing capacities */
|
|
cap_in; /* sum of incoming capacities */
|
|
|
|
long no_lines=0, /* no of current input line */
|
|
/* no_plines=0, */ /* no of problem-lines */
|
|
/* no_nlines=0, */ /* no of node lines */
|
|
no_alines=0, /* no of arc-lines */
|
|
pos_current=0; /* 2*no_alines */
|
|
|
|
int /* k, */ /* temporary */
|
|
err_no; /* no of detected error */
|
|
|
|
/* -------------- error numbers & error messages ---------------- */
|
|
#define EN1 0
|
|
#define EN2 1
|
|
#define EN3 2
|
|
#define EN4 3
|
|
#define EN6 4
|
|
#define EN10 5
|
|
#define EN7 6
|
|
#define EN8 7
|
|
#define EN9 8
|
|
#define EN11 9
|
|
#define EN12 10
|
|
#define EN13 11
|
|
#define EN14 12
|
|
#define EN16 13
|
|
#define EN15 14
|
|
#define EN17 15
|
|
#define EN18 16
|
|
#define EN21 17
|
|
#define EN19 18
|
|
#define EN20 19
|
|
#define EN22 20
|
|
|
|
static char *err_message[] =
|
|
{
|
|
/* 0*/ "more than one problem line",
|
|
/* 1*/ "wrong number of parameters in the problem line",
|
|
/* 2*/ "it is not a Min-cost problem line",
|
|
/* 3*/ "bad value of a parameter in the problem line",
|
|
/* 4*/ "can't obtain enough memory to solve this problem",
|
|
/* 5*/ "",
|
|
/* 6*/ "can't read problem name",
|
|
/* 7*/ "problem description must be before node description",
|
|
/* 8*/ "wrong capacity bounds",
|
|
/* 9*/ "wrong number of parameters in the node line",
|
|
/*10*/ "wrong value of parameters in the node line",
|
|
/*11*/ "unbalanced problem",
|
|
/*12*/ "node descriptions must be before arc descriptions",
|
|
/*13*/ "too many arcs in the input",
|
|
/*14*/ "wrong number of parameters in the arc line",
|
|
/*15*/ "wrong value of parameters in the arc line",
|
|
/*16*/ "unknown line type in the input",
|
|
/*17*/ "read error",
|
|
/*18*/ "not enough arcs in the input",
|
|
/*19*/ "warning: capacities too big - excess overflow possible",
|
|
/*20*/ "can't read anything from the input file",
|
|
/*21*/ "warning: infinite capacity replaced by BIGGEST_FLOW"
|
|
};
|
|
/* --------------------------------------------------------------- */
|
|
|
|
|
|
|
|
/* set up */
|
|
nnodes=nNrow*nNcol+1; /* add one for ground node */
|
|
narcs=2*((nNrow+1)*nNcol+nNrow*(nNcol+1)); /* 2x for two directional arcs */
|
|
cumsupply=0;
|
|
nresidues=0;
|
|
|
|
/* get memory (formerly case 'p' in DIMACS file read) */
|
|
fprintf(sp2,"Setting up data structures for cs2 MCF solver\n");
|
|
n=nnodes;
|
|
m=narcs;
|
|
if ( n <= 0 || m <= 0 )
|
|
/*wrong value of no of arcs or nodes*/
|
|
{ err_no = EN4; goto error; }
|
|
|
|
/* allocating memory for 'nodes', 'arcs' and internal arrays */
|
|
nodes = (node*) CAlloc ( n+2, sizeof(node) );
|
|
arcs = (arc*) CAlloc ( 2*m+1, sizeof(arc) );
|
|
cap = (short*) CAlloc ( 2*m, sizeof(short) ); /* changed to short */
|
|
arc_tail = (long*) CAlloc ( 2*m, sizeof(long) );
|
|
arc_first= (long*) CAlloc ( n+2, sizeof(long) );
|
|
/* arc_first [ 0 .. n+1 ] = 0 - initialized by calloc */
|
|
|
|
for ( in = nodes; in <= nodes + n; in ++ )
|
|
in -> excess = 0;
|
|
|
|
if ( nodes == NULL || arcs == NULL ||
|
|
arc_first == NULL || arc_tail == NULL )
|
|
/* memory is not allocated */
|
|
{ err_no = EN6; goto error; }
|
|
|
|
/* setting pointer to the first arc */
|
|
arc_current = arcs;
|
|
node_max = 0;
|
|
node_min = n;
|
|
m_c = 0;
|
|
total_p = total_n = 0;
|
|
|
|
for ( ndp = nodes; ndp < nodes + n; ndp ++ )
|
|
ndp -> excess = 0;
|
|
|
|
/* end of former case 'p' */
|
|
|
|
|
|
/* load supply/demand info into arrays (case 'n' in former loop) */
|
|
for(col=0; col<nNcol; col++){
|
|
for(row=0; row<nNrow; row++){
|
|
if(residue[row][col]){
|
|
i=(col*nNrow + row + 1);
|
|
excess=residue[row][col];
|
|
( nodes + i ) -> excess = excess;
|
|
if ( excess > 0 ) total_p += (double)excess;
|
|
if ( excess < 0 ) total_n -= (double)excess;
|
|
nresidues++;
|
|
cumsupply+=residue[row][col];
|
|
}
|
|
}
|
|
}
|
|
|
|
/* give ground node excess of -cumsupply */
|
|
( nodes + nnodes ) -> excess = -cumsupply;
|
|
if (cumsupply < 0) total_p -= (double)cumsupply;
|
|
if (cumsupply > 0) total_n += (double)cumsupply;
|
|
|
|
/* load arc info into arrays (case 'a' in former loop) */
|
|
low=0;
|
|
acap=ARCUBOUND;
|
|
|
|
/* horizontal (row) direction arcs first */
|
|
for(arcctr=1;arcctr<=2*nNrow*nNcol+nNrow+nNcol;arcctr++){
|
|
if(arcctr<=nNrow*(nNcol+1)){
|
|
/* row (horizontal) arcs first */
|
|
nodectr=arcctr;
|
|
if(nodectr<=nNrow*nNcol){
|
|
tail=nodectr;
|
|
}else{
|
|
tail=nnodes;
|
|
}
|
|
if(nodectr<=nNrow){
|
|
head=nnodes;
|
|
}else{
|
|
head=nodectr-nNrow;
|
|
}
|
|
cost=rowcost[((nodectr-1) % nNrow)][(int )((nodectr-1)/nNrow)];
|
|
}else{
|
|
/* column (vertical) arcs */
|
|
nodectr=arcctr-nNrow*(nNcol+1);
|
|
if(nodectr % (nNrow+1)==0){
|
|
tail=nnodes;
|
|
}else{
|
|
tail=(int )(nodectr-ceil(nodectr/(nNrow+1.0))+1);
|
|
}
|
|
if(nodectr % (nNrow+1)==1){
|
|
head=nnodes;
|
|
}else{
|
|
head=(int )(nodectr-ceil(nodectr/(nNrow+1.0)));
|
|
}
|
|
cost=colcost[((nodectr-1) % (nNrow+1))][(int )((nodectr-1)/(nNrow+1))];
|
|
}
|
|
|
|
if ( tail < 0 || tail > n ||
|
|
head < 0 || head > n
|
|
)
|
|
/* wrong value of nodes */
|
|
{ err_no = EN17; goto error; }
|
|
|
|
if ( acap < 0 ) {
|
|
acap = BIGGEST_FLOW;
|
|
if (!inf_cap) {
|
|
inf_cap = 1;
|
|
fprintf ( sp0, "\ncs2 solver: %s\n", err_message[21] );
|
|
}
|
|
}
|
|
|
|
if ( low < 0 || low > acap )
|
|
{ err_no = EN9; goto error; }
|
|
|
|
for(dir=0;dir<=1;dir++){
|
|
if(dir){
|
|
/* switch head and tail and loop for two directional arcs */
|
|
temp=tail;
|
|
tail=head;
|
|
head=temp;
|
|
}
|
|
|
|
/* no of arcs incident to node i is placed in arc_first[i+1] */
|
|
arc_first[tail + 1] ++;
|
|
arc_first[head + 1] ++;
|
|
in = nodes + tail;
|
|
jn = nodes + head;
|
|
dcost = (double)cost;
|
|
|
|
/* storing information about the arc */
|
|
arc_tail[pos_current] = tail;
|
|
arc_tail[pos_current+1] = head;
|
|
arc_current -> head = jn;
|
|
arc_current -> r_cap = acap - low;
|
|
cap[pos_current] = acap;
|
|
arc_current -> cost = dcost;
|
|
arc_current -> sister = arc_current + 1;
|
|
( arc_current + 1 ) -> head = nodes + tail;
|
|
( arc_current + 1 ) -> r_cap = 0;
|
|
cap[pos_current+1] = 0;
|
|
( arc_current + 1 ) -> cost = -dcost;
|
|
( arc_current + 1 ) -> sister = arc_current;
|
|
|
|
in -> excess -= low;
|
|
jn -> excess += low;
|
|
|
|
/* searching for minimum and maximum node */
|
|
if ( head < node_min ) node_min = head;
|
|
if ( tail < node_min ) node_min = tail;
|
|
if ( head > node_max ) node_max = head;
|
|
if ( tail > node_max ) node_max = tail;
|
|
|
|
if ( dcost < 0 ) dcost = -dcost;
|
|
if ( dcost > m_c && acap > 0 ) m_c = dcost;
|
|
|
|
no_alines ++;
|
|
arc_current += 2;
|
|
pos_current += 2;
|
|
|
|
}/* end of for loop over arc direction */
|
|
}/* end of for loop over arcss */
|
|
|
|
|
|
/* ----- all is red or error while reading ----- */
|
|
|
|
if ( ABS( total_p - total_n ) > 0.5 ) /* unbalanced problem */
|
|
{ err_no = EN13; goto error; }
|
|
|
|
/********** ordering arcs - linear time algorithm ***********/
|
|
|
|
/* first arc from the first node */
|
|
( nodes + node_min ) -> first = arcs;
|
|
|
|
/* before below loop arc_first[i+1] is the number of arcs outgoing from i;
|
|
after this loop arc_first[i] is the position of the first
|
|
outgoing from node i arcs after they would be ordered;
|
|
this value is transformed to pointer and written to node.first[i]
|
|
*/
|
|
|
|
for ( i = node_min + 1; i <= node_max + 1; i ++ )
|
|
{
|
|
arc_first[i] += arc_first[i-1];
|
|
( nodes + i ) -> first = arcs + arc_first[i];
|
|
}
|
|
|
|
|
|
for ( i = node_min; i < node_max; i ++ ) /* scanning all the nodes
|
|
exept the last*/
|
|
{
|
|
|
|
last = ( ( nodes + i + 1 ) -> first ) - arcs;
|
|
/* arcs outgoing from i must be cited
|
|
from position arc_first[i] to the position
|
|
equal to initial value of arc_first[i+1]-1 */
|
|
|
|
for ( arc_num = arc_first[i]; arc_num < last; arc_num ++ )
|
|
{ tail = arc_tail[arc_num];
|
|
|
|
while ( tail != i )
|
|
/* the arc no arc_num is not in place because arc cited here
|
|
must go out from i;
|
|
we'll put it to its place and continue this process
|
|
until an arc in this position would go out from i */
|
|
|
|
{ arc_new_num = arc_first[tail];
|
|
arc_current = arcs + arc_num;
|
|
arc_new = arcs + arc_new_num;
|
|
|
|
/* arc_current must be cited in the position arc_new
|
|
swapping these arcs: */
|
|
|
|
head_p = arc_new -> head;
|
|
arc_new -> head = arc_current -> head;
|
|
arc_current -> head = head_p;
|
|
|
|
acap = cap[arc_new_num];
|
|
cap[arc_new_num] = cap[arc_num];
|
|
cap[arc_num] = acap;
|
|
|
|
acap = arc_new -> r_cap;
|
|
arc_new -> r_cap = arc_current -> r_cap;
|
|
arc_current -> r_cap = acap;
|
|
|
|
dcost = arc_new -> cost;
|
|
arc_new -> cost = arc_current -> cost;
|
|
arc_current -> cost = dcost;
|
|
|
|
if ( arc_new != arc_current -> sister )
|
|
{
|
|
arc_tmp = arc_new -> sister;
|
|
arc_new -> sister = arc_current -> sister;
|
|
arc_current -> sister = arc_tmp;
|
|
|
|
( arc_current -> sister ) -> sister = arc_current;
|
|
( arc_new -> sister ) -> sister = arc_new;
|
|
}
|
|
|
|
arc_tail[arc_num] = arc_tail[arc_new_num];
|
|
arc_tail[arc_new_num] = tail;
|
|
|
|
/* we increase arc_first[tail] */
|
|
arc_first[tail] ++ ;
|
|
|
|
tail = arc_tail[arc_num];
|
|
}
|
|
}
|
|
/* all arcs outgoing from i are in place */
|
|
}
|
|
|
|
/* ----------------------- arcs are ordered ------------------------- */
|
|
|
|
/*------------ testing network for possible excess overflow ---------*/
|
|
|
|
for ( ndp = nodes + node_min; ndp <= nodes + node_max; ndp ++ )
|
|
{
|
|
cap_in = ( ndp -> excess );
|
|
cap_out = - ( ndp -> excess );
|
|
for ( arc_current = ndp -> first; arc_current != (ndp+1) -> first;
|
|
arc_current ++ )
|
|
{
|
|
arc_num = arc_current - arcs;
|
|
if ( cap[arc_num] > 0 ) cap_out += cap[arc_num];
|
|
if ( cap[arc_num] == 0 )
|
|
cap_in += cap[( arc_current -> sister )-arcs];
|
|
}
|
|
|
|
/*
|
|
if (cap_in > BIGGEST_FLOW || cap_out > BIGGEST_FLOW)
|
|
{
|
|
fprintf ( sp0, "\ncs2 solver: %s\n", err_message[EN20] );
|
|
break;
|
|
}
|
|
*/
|
|
}
|
|
|
|
/* ----------- assigning output values ------------*/
|
|
*m_ad = m;
|
|
*n_ad = node_max - node_min + 1;
|
|
*node_min_ad = node_min;
|
|
*nodes_ad = nodes + node_min;
|
|
*arcs_ad = arcs;
|
|
*m_c_ad = m_c;
|
|
*cap_ad = cap;
|
|
|
|
/* free internal memory */
|
|
free ( arc_first ); free ( arc_tail );
|
|
|
|
/* Thanks God! All is done! */
|
|
return (0);
|
|
|
|
/* ---------------------------------- */
|
|
error: /* error found reading input */
|
|
|
|
fprintf ( sp0, "\ncs2 solver: line %ld of input - %s\n",
|
|
no_lines, err_message[err_no] );
|
|
|
|
exit (ABNORMAL_EXIT);
|
|
|
|
/* this is a needless return statement so the compiler doesn't complain */
|
|
return(1);
|
|
|
|
}
|
|
/* -------------------- end of parser -------------------*/
|
|
|
|
|
|
|