Thread View: comp.sources.reviewed
1 messages
1 total messages
Started by Michel Berkelaar
Wed, 05 Aug 1992 22:50
v02i023: lp_solve - linear problem solver, Part02/02
Author: Michel Berkelaar
Date: Wed, 05 Aug 1992 22:50
Date: Wed, 05 Aug 1992 22:50
1672 lines
48011 bytes
48011 bytes
Submitted-by: Michel Berkelaar <michel@es.ele.tue.nl>
Posting-number: Volume 2, Issue 23
Archive-name: lp_solve/part02
#! /bin/sh
# This is a shell archive. Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file". To overwrite existing
# files, type "sh file -c". You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g.. If this archive is complete, you
# will see the following message at the end:
# "End of archive 2 (of 2)."
# Contents: lp_solve/ex5.lp lp_solve/read.c lp_solve/solve.c
# Wrapped by michel@varda.es.ele.tue.nl on Fri Jun 19 09:05:17 1992
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f lp_solve/ex5.lp -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"lp_solve/ex5.lp\"
else
echo shar: Extracting \"lp_solve/ex5.lp\" \(8873 characters\)
sed "s/^X//" >lp_solve/ex5.lp <<'END_OF_lp_solve/ex5.lp'
X/* correct value of objective function after optimization: 2.3314e+06 */
X
X + 35 COL00001 - 54.9 COL00002 - 54.9 COL00003 + 35 COL00004
X - 23.5 COL00005 - 23.5 COL00006 - 8.72 COL00007 - 9.72 COL00008
X - 6.74 COL00009 - 6.84 COL00010 - 15 COL00011 - 22.5 COL00012 + 500 COL00014
X - 18.7 COL00015 + 258.3 COL00016 + 662 COL00017 - 3 COL00018 - 0.39 COL00019
X - 0.47 COL00020 + 35 COL00021 - 54.9 COL00022 - 54.9 COL00023 + 35 COL00024
X - 23.5 COL00025 - 23.5 COL00026 - 8.72 COL00027 - 9.72 COL00028
X - 6.74 COL00029 - 6.84 COL00030 - 15 COL00031 - 22.5 COL00032 + 500 COL00034
X - 18.7 COL00035 + 258.3 COL00036 + 662 COL00037 - 3 COL00038 - 0.39 COL00039
X - 0.47 COL00040 + 35 COL00041 - 54.9 COL00042 - 54.9 COL00043 + 35 COL00044
X - 23.5 COL00045 - 23.5 COL00046 - 8.72 COL00047 - 9.72 COL00048
X - 6.74 COL00049 - 6.84 COL00050 - 15 COL00051 - 22.5 COL00052 + 500 COL00054
X - 18.7 COL00055 + 258.3 COL00056 + 662 COL00057 - 3 COL00058 - 0.39 COL00059
X - 0.47 COL00060 + 35 COL00061 - 54.9 COL00062 - 54.9 COL00063 + 35 COL00064
X - 23.5 COL00065 - 23.5 COL00066 - 8.72 COL00067 - 9.72 COL00068
X - 6.74 COL00069 - 6.84 COL00070 - 15 COL00071 - 22.5 COL00072 + 500 COL00074
X - 18.7 COL00075 + 258.3 COL00076 + 662 COL00077 - 3 COL00078 - 0.39 COL00079
X - 0.47 COL00080 + 35 COL00081 - 54.9 COL00082 - 54.9 COL00083 + 35 COL00084
X - 23.5 COL00085 - 23.5 COL00086 - 8.72 COL00087 - 9.72 COL00088
X - 6.74 COL00089 - 6.84 COL00090 - 15 COL00091 - 22.5 COL00092 + 500 COL00094
X - 18.7 COL00095 + 258.3 COL00096 + 662 COL00097 - 3 COL00098 - 0.39 COL00099
X - 0.47 COL00100 + 35 COL00101 - 54.9 COL00102 - 54.9 COL00103 + 35 COL00104
X - 23.5 COL00105 - 23.5 COL00106 - 8.72 COL00107 - 9.72 COL00108
X - 6.74 COL00109 - 6.84 COL00110 - 15 COL00111 - 22.5 COL00112 + 500 COL00114
X - 18.7 COL00115 + 258.3 COL00116 + 662 COL00117 - 3 COL00118 - 0.39 COL00119
X - 0.47 COL00120 + 35 COL00121 - 54.9 COL00122 - 54.9 COL00123 + 35 COL00124
X - 23.5 COL00125 - 23.5 COL00126 - 8.72 COL00127 - 9.72 COL00128
X - 6.74 COL00129 - 6.84 COL00130 - 15 COL00131 - 22.5 COL00132 + 500 COL00134
X - 18.7 COL00135 + 258.3 COL00136 + 662 COL00137 - 3 COL00138 - 0.39 COL00139
X - 0.47 COL00140 ;
X COL00001 + COL00002 = 158 ;
X COL00003 + COL00004 = 158 ;
X COL00002 - COL00005 = 0 ;
X COL00003 - COL00006 = 0 ;
X - 0.4 COL00002 - 0.4 COL00003 - 1.7 COL00005 - 1.7 COL00006
X + 1.5 COL00008 + COL00011 = 3092.96 ;
X COL00008 <= 2566.67 ;
X COL00007 + COL00008 + COL00009 + COL00010 + COL00013 = 6900 ;
X COL00013 <= 1600 ;
X COL00009 + COL00010 >= 800 ;
X COL00005 - COL00014 - COL00015 = 0 ;
X COL00006 - COL00016 = 0 ;
X - COL00017 = -375.2 ;
X - 0.245 COL00005 - 0.245 COL00006 + 0.2 COL00007 - COL00018
X = -92.12 ;
X 1.2 COL00009 - COL00019 = -684 ;
X - 0.5 COL00002 - 0.5 COL00003 - 1.4 COL00005 - 1.4 COL00006
X + COL00010 + COL00012 - COL00020 = -150 ;
X - 0.5 COL00017 + COL00021 + COL00022 = 0 ;
X - 0.5 COL00017 + COL00023 + COL00024 = 0 ;
X COL00022 - COL00025 = 0 ;
X COL00023 - COL00026 = 0 ;
X COL00018 <= 1800 ;
X COL00019 <= 2400 ;
X COL00020 <= 1200 ;
X - 0.48 COL00014 - 0.48 COL00015 - 0.48 COL00016 - 9.32 COL00017
X - 0.4 COL00022 - 0.4 COL00023 - 1.7 COL00025 - 1.7 COL00026
X + 1.5 COL00028 + COL00031 = 0 ;
X COL00028 <= 2566.67 ;
X COL00027 + COL00028 + COL00029 + COL00030 + COL00033 = 6900 ;
X COL00033 <= 1600 ;
X COL00029 + COL00030 >= 800 ;
X COL00025 - COL00034 - COL00035 = 0 ;
X COL00026 - COL00036 = 0 ;
X COL00015 + 0.7 COL00017 - COL00037 = 0 ;
X - 0.49 COL00014 - 0.49 COL00015 - 0.49 COL00016 - 0.56 COL00017
X + COL00018 - 0.245 COL00025 - 0.245 COL00026 + 0.2 COL00027
X - COL00038 = 0 ;
X - COL00017 + COL00019 + 1.2 COL00029 - COL00039 = 0 ;
X COL00020 - 0.5 COL00022 - 0.5 COL00023 - 1.4 COL00025 - 1.4 COL00026
X + COL00030 + COL00032 - COL00040 = 0 ;
X COL00017 - COL00037 <= 0 ;
X - 0.5 COL00037 + COL00041 + COL00042 = 0 ;
X - 0.5 COL00037 + COL00043 + COL00044 = 0 ;
X COL00042 - COL00045 = 0 ;
X COL00043 - COL00046 = 0 ;
X COL00038 <= 1800 ;
X COL00039 <= 2400 ;
X COL00040 <= 1200 ;
X - 0.48 COL00034 - 0.48 COL00035 - 0.48 COL00036 - 9.32 COL00037
X - 0.4 COL00042 - 0.4 COL00043 - 1.7 COL00045 - 1.7 COL00046
X + 1.5 COL00048 + COL00051 = 0 ;
X COL00048 <= 2566.67 ;
X COL00047 + COL00048 + COL00049 + COL00050 + COL00053 = 6900 ;
X COL00053 <= 1600 ;
X COL00049 + COL00050 >= 800 ;
X COL00045 - COL00054 - COL00055 = 0 ;
X COL00046 - COL00056 = 0 ;
X COL00035 + 0.7 COL00037 - COL00057 = 0 ;
X - 0.49 COL00034 - 0.49 COL00035 - 0.49 COL00036 - 0.56 COL00037
X + COL00038 - 0.245 COL00045 - 0.245 COL00046 + 0.2 COL00047
X - COL00058 = 0 ;
X - COL00037 + COL00039 + 1.2 COL00049 - COL00059 = 0 ;
X COL00040 - 0.5 COL00042 - 0.5 COL00043 - 1.4 COL00045 - 1.4 COL00046
X + COL00050 + COL00052 - COL00060 = 0 ;
X COL00037 - COL00057 <= 0 ;
X - 0.5 COL00057 + COL00061 + COL00062 = 0 ;
X - 0.5 COL00057 + COL00063 + COL00064 = 0 ;
X COL00062 - COL00065 = 0 ;
X COL00063 - COL00066 = 0 ;
X COL00058 <= 1800 ;
X COL00059 <= 2400 ;
X COL00060 <= 1200 ;
X - 0.48 COL00054 - 0.48 COL00055 - 0.48 COL00056 - 9.32 COL00057
X - 0.4 COL00062 - 0.4 COL00063 - 1.7 COL00065 - 1.7 COL00066
X + 1.5 COL00068 + COL00071 = 0 ;
X COL00068 <= 2566.67 ;
X COL00067 + COL00068 + COL00069 + COL00070 + COL00073 = 6900 ;
X COL00073 <= 1600 ;
X COL00069 + COL00070 >= 800 ;
X COL00065 - COL00074 - COL00075 = 0 ;
X COL00066 - COL00076 = 0 ;
X COL00055 + 0.7 COL00057 - COL00077 = 0 ;
X - 0.49 COL00054 - 0.49 COL00055 - 0.49 COL00056 - 0.56 COL00057
X + COL00058 - 0.245 COL00065 - 0.245 COL00066 + 0.2 COL00067
X - COL00078 = 0 ;
X - COL00057 + COL00059 + 1.2 COL00069 - COL00079 = 0 ;
X COL00060 - 0.5 COL00062 - 0.5 COL00063 - 1.4 COL00065 - 1.4 COL00066
X + COL00070 + COL00072 - COL00080 = 0 ;
X COL00057 - COL00077 <= 0 ;
X - 0.5 COL00077 + COL00081 + COL00082 = 0 ;
X - 0.5 COL00077 + COL00083 + COL00084 = 0 ;
X COL00082 - COL00085 = 0 ;
X COL00083 - COL00086 = 0 ;
X COL00078 <= 1800 ;
X COL00079 <= 2400 ;
X COL00080 <= 1200 ;
X - 0.48 COL00074 - 0.48 COL00075 - 0.48 COL00076 - 9.32 COL00077
X - 0.4 COL00082 - 0.4 COL00083 - 1.7 COL00085 - 1.7 COL00086
X + 1.5 COL00088 + COL00091 = 0 ;
X COL00088 <= 2566.67 ;
X COL00087 + COL00088 + COL00089 + COL00090 + COL00093 = 6900 ;
X COL00093 <= 1600 ;
X COL00089 + COL00090 >= 800 ;
X COL00085 - COL00094 - COL00095 = 0 ;
X COL00086 - COL00096 = 0 ;
X COL00075 + 0.7 COL00077 - COL00097 = 0 ;
X - 0.49 COL00074 - 0.49 COL00075 - 0.49 COL00076 - 0.56 COL00077
X + COL00078 - 0.245 COL00085 - 0.245 COL00086 + 0.2 COL00087
X - COL00098 = 0 ;
X - COL00077 + COL00079 + 1.2 COL00089 - COL00099 = 0 ;
X COL00080 - 0.5 COL00082 - 0.5 COL00083 - 1.4 COL00085 - 1.4 COL00086
X + COL00090 + COL00092 - COL00100 = 0 ;
X COL00077 - COL00097 <= 0 ;
X - 0.5 COL00097 + COL00101 + COL00102 = 0 ;
X - 0.5 COL00097 + COL00103 + COL00104 = 0 ;
X COL00102 - COL00105 = 0 ;
X COL00103 - COL00106 = 0 ;
X COL00098 <= 1800 ;
X COL00099 <= 2400 ;
X COL00100 <= 1200 ;
X - 0.48 COL00094 - 0.48 COL00095 - 0.48 COL00096 - 9.32 COL00097
X - 0.4 COL00102 - 0.4 COL00103 - 1.7 COL00105 - 1.7 COL00106
X + 1.5 COL00108 + COL00111 = 0 ;
X COL00108 <= 2566.67 ;
X COL00107 + COL00108 + COL00109 + COL00110 + COL00113 = 6900 ;
X COL00113 <= 1600 ;
X COL00109 + COL00110 >= 800 ;
X COL00105 - COL00114 - COL00115 = 0 ;
X COL00106 - COL00116 = 0 ;
X COL00095 + 0.7 COL00097 - COL00117 = 0 ;
X - 0.49 COL00094 - 0.49 COL00095 - 0.49 COL00096 - 0.56 COL00097
X + COL00098 - 0.245 COL00105 - 0.245 COL00106 + 0.2 COL00107
X - COL00118 = 0 ;
X - COL00097 + COL00099 + 1.2 COL00109 - COL00119 = 0 ;
X COL00100 - 0.5 COL00102 - 0.5 COL00103 - 1.4 COL00105 - 1.4 COL00106
X + COL00110 + COL00112 - COL00120 = 0 ;
X COL00097 - COL00117 <= 0 ;
X - 0.5 COL00117 + COL00121 + COL00122 = 0 ;
X - 0.5 COL00117 + COL00123 + COL00124 = 0 ;
X COL00122 - COL00125 = 0 ;
X COL00123 - COL00126 = 0 ;
X COL00118 <= 1800 ;
X COL00119 <= 2400 ;
X COL00120 <= 1200 ;
X - 0.48 COL00114 - 0.48 COL00115 - 0.48 COL00116 - 9.32 COL00117
X - 0.4 COL00122 - 0.4 COL00123 - 1.7 COL00125 - 1.7 COL00126
X + 1.5 COL00128 + COL00131 = 0 ;
X COL00128 <= 2566.67 ;
X COL00127 + COL00128 + COL00129 + COL00130 + COL00133 = 6900 ;
X COL00133 <= 1600 ;
X COL00129 + COL00130 >= 800 ;
X COL00125 - COL00134 - COL00135 = 0 ;
X COL00126 - COL00136 = 0 ;
X COL00115 + 0.7 COL00117 - COL00137 = 0 ;
X - 0.49 COL00114 - 0.49 COL00115 - 0.49 COL00116 - 0.56 COL00117
X + COL00118 - 0.245 COL00125 - 0.245 COL00126 + 0.2 COL00127
X - COL00138 = 0 ;
X - COL00117 + COL00119 + 1.2 COL00129 - COL00139 = 0 ;
X COL00120 - 0.5 COL00122 - 0.5 COL00123 - 1.4 COL00125 - 1.4 COL00126
X + COL00130 + COL00132 - COL00140 = 0 ;
X COL00117 - COL00137 <= 0 ;
X
END_OF_lp_solve/ex5.lp
if test 8873 -ne `wc -c <lp_solve/ex5.lp`; then
echo shar: \"lp_solve/ex5.lp\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f lp_solve/read.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"lp_solve/read.c\"
else
echo shar: Extracting \"lp_solve/read.c\" \(13928 characters\)
sed "s/^X//" >lp_solve/read.c <<'END_OF_lp_solve/read.c'
X/*
X============================================================================
XNAME : read.c
XPURPOSE : translation of lp-problem and storage in sparse matrix
XSHORT : Subroutines for yacc program to store the input in an intermediate
X data-structure. The yacc and lex programs translate the input.
X First the problemsize is determined and the date is read into
X an intermediate structure, then readinput fills the sparse matrix.
XUSAGE : call yyparse(); to start reading the input.
X call readinput(); to fill the sparse matrix.
X============================================================================
XRows : contains the amount of rows + 1
X Rows-1 is the amount of constraints (no bounds)
X Rows also contains the rownr 0 which is the objectfunction
XColumns : contains the amount of columns (different variable names
X found in the constraints)
XNonnuls : contains the amount of nonnuls = sum of different entries
X of all columns in the constraints and in the objectfunction
XHash_tab : contains all columnnames on the first level of the structure
X the row information is kept under each column structure
X in a linked list (also the objext funtion is in this structure)
X Bound information is also stored under under the column name
XFirst_rside : points to a linked list containing all relational operators
X and the righthandside values of the constraints
X the linked list is in reversed order with respect to the
X rownumbers
X============================================================================
X*/
X#include "defines.h"
X#include "globals.h"
X
X/*---------------------------------------------------------------------------*/
X
X/*
X * errorhandeling routine for yyparse()
X */
Xvoid yyerror(char *string)
X{
X fprintf(stderr, "PARSING ERROR: %s on line %d, quiting\n", string, yylineno);
X exit(1);
X}
X
Xvoid check_decl(char *str)
X{
X if(strcmp(str, "int"))
X {
X fprintf(stderr, "Unknown declaration specifier %s on line %d, ignored\n",
X str, yylineno);
X Ignore_decl = TRUE;
X }
X}
X
Xvoid add_int_var(char *name)
X{
X hashelem *hp;
X
X if(Verbose)
X fprintf(stderr, "int: %s\n", name);
X if(!(hp = gethash(name)))
X fprintf(stderr,
X "Unknown variable %s declared integer on line %d, ignored\n",
X name, yylineno);
X else
X hp->must_be_int = 1;
X}
X
X
X/*
X * initialisation of hashtable and globals.
X */
Xvoid init_read(void)
X{
X int i;
X Rows = 0;
X Nonnuls = 0;
X Columns = 0;
X for(i = 0; i<HASH_SIZE; Hash_tab[i++] = NULL);
X CALLOC(First_rside, 1, rside);
X First_rside->value = 0;
X First_rside->relat = OF; /* first row (nr 0) is always the objective function */
X} /* init */
X
X/*---------------------------------------------------------------------------*/
X
X/*
X * Returns hashvalue of string
X * hashvalue = sum ( asciivalue of character * indexnumber) mod HASHSIZE
X */
Xint hashvalue(char *string)
X{
X int i, j;
X i = j = 0;
X while((int)string[j]&&j<MAXSTRL)
X i += (NAMELEN-j)*(int)string[j++];
X return(i % HASH_SIZE);
X} /* hashvalue */
X
X
X
X/*
X * Returns a pointer to hashelement with colname = variable
X * If this hashelement does not exists, gethash() returns a NULL pointer
X */
Xhashelem *gethash(char *variable)
X{
X hashelem *h_tab_p;
X for(h_tab_p = Hash_tab[hashvalue(variable)];
X h_tab_p != NULL;
X h_tab_p = h_tab_p->next)
X if(strcmp(variable, h_tab_p->colname) == 0)
X return(h_tab_p);
X return(h_tab_p);
X} /* gethash */
X
X
X
X/*
X * searchs in column-list (p is pointer to first element of column-list)
X * for column->row = row.
X * getrow() returns a pointer to this column structure.
X * If not found a NULL-pointer is returned
X */
Xcolumn *getrow(column *p,
X int row)
X{
X for(; p != NULL; p = p->next)
X if(p->row == row)
X return(p);
X return(p) ;
X} /* getrow */
X
X
X
X/*
X * Creates a bound record.
X * Set lowbo = 0 and upbo = INFINITE
X *
X */
Xbound *create_bound_rec(void)
X{
X bound *bp;
X CALLOC(bp, 1, bound);
X bp->upbo = INFINITE;
X bp->lowbo = 0;
X return(bp);
X} /* create_bound_rec */
X
X
X
X/*
X * clears the tmp_store variable after all information has been copied
X */
Xvoid null_tmp_store(void)
X{
X tmp_store.value = 0;
X tmp_store.rhs_value = 0;
X}
X
X/*---------------------------------------------------------------------------*/
X
X/*
X * variable : pointer to text array with name of variable
X * row : the rownumber of the constraint
X * value : value of matrixelement
X * A(row, variable).
X * Sign : (global) determines the sign of value.
X * store() : stores value in matrix
X * A(row, variable). If A(row, variable) already contains data,
X * value is added to the existing value.
X */
Xvoid store(char *variable,
X int row,
X double value)
X{
X hashelem *h_tab_p;
X column *col_p;
X int hv;
X
X if(value == (double)0)
X return;
X if((h_tab_p = gethash(variable)) == NULL)
X {
X CALLOC(h_tab_p, 1, hashelem);
X Columns++; /* counter for calloc of final array */
X hv = hashvalue(variable);
X h_tab_p->next = Hash_tab[hv];
X Hash_tab[hv] = h_tab_p;
X strcpy(h_tab_p->colname, variable);
X CALLOC(h_tab_p->col, 1, column);
X Nonnuls++; /* for calloc of final arrays */
X h_tab_p->col->row = row;
X h_tab_p->col->value = value;
X }
X else
X if((col_p = getrow(h_tab_p->col, row)) == NULL)
X {
X CALLOC(col_p, 1, column);
X Nonnuls++; /* for calloc of final arrays */
X col_p->value = value;
X col_p->row = row;
X col_p->next = h_tab_p->col;
X h_tab_p->col = col_p;
X }
X else
X col_p->value += value;
X} /* store */
X
X
X
X/*---------------------------------------------------------------------------*/
X
X/*
X * store relational operator given in yylex[0] in the rightside list.
X * Also checks if it constaint was a bound and if so stores it in the
X * boundslist
X */
Xvoid store_re_op(void)
X{
X short tmp_relat;
X switch(yytext[0]) {
X
X case '=':
X tmp_relat = EQ;
X break;
X
X case '>':
X tmp_relat=GE;
X break;
X
X case '<':
X tmp_relat=LE;
X break;
X
X default:
X break;
X }
X
X if(Lin_term_count > 1) /* it is not a bound */
X First_rside->relat = tmp_relat;
X else /* could be a bound */
X tmp_store.relat = tmp_relat;
X} /* save_re_op */
X
X
X
X/*
X * store RHS value in the rightside structure
X * if type = true then
X */
Xvoid rhs_store(double value)
X{
X if(Lin_term_count <= 1) /* could be a bound */
X tmp_store.rhs_value += value;
X else /* not a bound */
X First_rside->value += value;
X} /* RHS_store */
X
X
X
X/*
X * store all data in the right place
X * count the amount of lineair terms in a constraint
X * only store in data-structure if the constraint is not a bound
X */
Xvoid var_store(char *var,
X int row,
X double value)
X{
X if(strlen(var) > MAXSTRL)
X {
X fprintf(stderr, "Variable name too long, at most %d characters allowed",
X MAXSTRL);
X exit(1);
X }
X /* also in a bound the same var name can occur more than once. Check for
X this. Don't increment Lin_term_count */
X
X if(Lin_term_count != 1 || strcmp(tmp_store.name, var) != 0)
X Lin_term_count++;
X
X if(row == 0) /* always store objective function with rownr == 0 */
X {
X store(var, row, value);
X return;
X }
X
X if(Lin_term_count == 1) /* don't yet store. could be a bound */
X {
X strcpy(tmp_store.name, var);
X tmp_store.row = row;
X tmp_store.value += value;
X return;
X }
X
X if(Lin_term_count == 2) /* now you can also store the first variable */
X {
X rside *rp;
X /* make space for the rhs information */
X CALLOC(rp, 1, rside);
X rp->next = First_rside;
X First_rside = rp;
X First_rside->value = tmp_store.rhs_value;
X First_rside->relat = tmp_store.relat;
X
X if (tmp_store.value != 0)
X store(tmp_store.name, tmp_store.row, tmp_store.value);
X
X null_tmp_store();
X }
X
X store(var, row, value);
X} /* var_store */
X
X
X
X/*
X * store the information in tmp_store because it is a bound
X */
Xstore_bounds(void)
X{
X if (tmp_store.value != 0)
X {
X hashelem *h_tab_p;
X int hv;
X
X if((h_tab_p = gethash(tmp_store.name)) == NULL)
X {
X /* a new columnname is found, create an entry in the hashlist */
X CALLOC(h_tab_p, 1, hashelem);
X Columns++; /* counter for calloc of final array */
X hv = hashvalue(tmp_store.name);
X h_tab_p->next = Hash_tab[hv];
X Hash_tab[hv] = h_tab_p;
X strcpy(h_tab_p->colname, tmp_store.name);
X /* create a place to store bounds information */
X h_tab_p->bnd = create_bound_rec();
X }
X else
X if(h_tab_p->bnd == NULL)
X /* create a place to store bounds information */
X h_tab_p->bnd = create_bound_rec();
X /* else bound_rec already exists */
X
X if(tmp_store.value < 0) /* divide by negative number, */
X /* relational operator may change */
X {
X if (tmp_store.relat == GE)
X tmp_store.relat = LE;
X else if (tmp_store.relat == LE)
X tmp_store.relat = GE;
X }
X
X if ((tmp_store.relat == GE) || (tmp_store.relat == EQ))
X h_tab_p->bnd->lowbo = tmp_store.rhs_value / tmp_store.value;
X if ((tmp_store.relat == LE) || (tmp_store.relat == EQ))
X h_tab_p->bnd->upbo = tmp_store.rhs_value / tmp_store.value;
X }
X
X null_tmp_store();
X} /* store_bounds */
X
X/* ========================================================================= */
X
X
X/*
X * reallocates eta
X */
Xvoid resize_eta(void)
X{
X Cur_eta_size *= 2;
X if (Verbose)
X printf("Resizing Eta_value and Eta_rownr, new size is %d\n", Cur_eta_size);
X if(!(Eta_value = realloc(Eta_value, Cur_eta_size * sizeof(double))))
X {
X fprintf(stderr, "Error: cannot realloc Eta_value to size %d (entries)\n",
X Cur_eta_size);
X exit(1);
X }
X if(!(Eta_rownr = realloc(Eta_rownr, Cur_eta_size * sizeof(int))))
X {
X fprintf(stderr, "Error: cannot realloc Eta_rownr to size %d (entries)\n",
X Cur_eta_size);
X exit(1);
X }
X} /* resize_eta */
X
X
X
X
X/*
X * transport the data from the intermediate structure to the sparse matrix
X * and free the intermediate structure
X */
Xvoid readinput(int *cend,
X double *rh,
X short *relat,
X double *lowbo,
X double *upbo,
X matrec *mat,
X nstring *names)
X{
X int i, j, k, index, nn_ind;
X column *cp,*tcp; /* tcp (temporary cp) points to memory-space to free */
X hashelem *hp,*thp;
X bound *bp;
X int x;
X rside *rp;
X intrec *irp;
X
X /* initialize lower and upper bound arrays */
X for (i = 0; i <= Sum; i++)
X {
X lowbo[i] = 0;
X upbo[i] = INFINITE;
X }
X
X /* fill names with the rownames */
X for (i = 0; i <= Rows; i++)
X {
X /* the first row is row zero (the objective function) */
X sprintf(names[i], "%s%d", STD_ROW_NAME_PREFIX, i);
X }
X
X for (i = Rows;i >= 0;i--)
X {
X rp = First_rside;
X relat[i] = rp->relat;
X rh[i] = rp->value;
X First_rside = rp->next;
X free(rp); /* free memory when data has been read */
X }
X
X /* change upperbound to zero if the relational operator is the equal sign */
X for (i = 1; i <= Rows; i++)
X if (relat[i] == EQ)
X upbo[i] = 0;
X
X /* start reading the Hash_list structure */
X index = 0;
X nn_ind = 0;
X
X for (i = 0;i < HASH_SIZE;i++)
X {
X hp = Hash_tab[i];
X while (hp != NULL)
X {
X /* put an index in the cend array when a new name is found */
X cend[index++] = nn_ind;
X
X /* check if it must be an integer variable */
X if(hp->must_be_int)
X {
X CALLOC(irp, 1, intrec);
X irp->varnr = Rows + index;
X irp->next = First_int;
X First_int = irp;
X }
X /* check for bound */
X if (hp->bnd != NULL)
X {
X bp = hp->bnd;
X lowbo[Rows+index] = bp->lowbo;
X upbo[Rows+index] = bp->upbo;
X free(bp); /* free memory when data has been read*/
X }
X
X /* copy name of column variable */
X strcpy(names[Rows+index], hp->colname);
X
X /* put matrix values in sparse matrix */
X cp = hp->col;
X while (cp!=NULL)
X {
X mat[nn_ind].rownr = cp->row;
X mat[nn_ind].value = cp->value;
X nn_ind++;
X tcp = cp;
X cp = cp->next;
X free(tcp); /* free memory when data has been read */
X }
X thp = hp;
X hp = hp->next;
X free(thp); /* free memory when data has been read */
X }
X }
X cend[index] = nn_ind;
X
X if (Verbose)
X {
X printf("\n");
X printf("**********Data read**********\n");
X printf("Rows : %d\n", Rows);
X printf("Columns : %d\n", Columns);
X printf("Nonnuls : %d\n", Nonnuls);
X printf("\nSparse Matrix\nRow right hand side\n");
X printf("%4s %8s %3s %s\n\n", "nr", "row", "rel", "Value");
X for (i = 0; i <= Rows; i++)
X {
X printf("%4d %8s ", i, names[i]);
X if (relat[i] == LE) printf(" < ");
X if (relat[i] == EQ) printf(" = ");
X if (relat[i] == GE) printf(" > ");
X if (relat[i] == OF) printf("ObF");
X printf(" %f\n", rh[i]);
X }
X
X printf("\nMatrix contents\n%8s %8s %s\n\n","colname","row","value");
X j = 0;
X for (i = 0; i < Nonnuls; i++)
X {
X if (i == cend[j])
X printf("%8s %8s %f\n", names[Rows+ ++j], names[mat[i].rownr],
X mat[i].value);
X else
X printf(" %8s %e\n", names[mat[i].rownr], mat[i].value);
X }
X printf("\nBounds\n%8s %3s %s\n\n", "name", "rel", "value");
X for (i = 0; i <= Sum; i++)
X {
X if (upbo[i] < (INFINITE/10000)) /* double is to small to contain
X INFINITE thus is rounded */
X printf("%8s < %f\n", names[i], upbo[i]);
X if (lowbo[i] > 0)
X printf("%8s > %f\n", names[i], lowbo[i]);
X }
X
X printf("\n**********End data**********\n\n");
X }
X} /* readinput */
X
X
X/* ===================== END OF read.c ===================================== */
X
END_OF_lp_solve/read.c
if test 13928 -ne `wc -c <lp_solve/read.c`; then
echo shar: \"lp_solve/read.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f lp_solve/solve.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"lp_solve/solve.c\"
else
echo shar: Extracting \"lp_solve/solve.c\" \(18876 characters\)
sed "s/^X//" >lp_solve/solve.c <<'END_OF_lp_solve/solve.c'
X#include "defines.h"
X#include "globals.h"
X
X
Xvoid condensecol(int rownr,
X int numc,
X double *pcol)
X{
X int i, elnr;
X
X if (Verbose)
X printf("condensecol\n");
X elnr = Endetacol[numc];
X
X if (elnr + Rows + 2 > Cur_eta_size) /* maximum local growth of Eta */
X resize_eta();
X
X for (i = 0; i <= Rows; i++)
X if (i != rownr && pcol[i] != 0)
X {
X Eta_rownr[elnr] = i;
X Eta_value[elnr] = pcol[i];
X elnr++;
X }
X Eta_rownr[elnr] = rownr;
X Eta_value[elnr] = pcol[rownr];
X elnr++;
X Endetacol[numc + 1] = elnr;
X} /* condensecol */
X
Xvoid addetacol(int *numc)
X{
X int i, j, k;
X double theta;
X
X if (Verbose)
X printf("addetacol\n");
X j = Endetacol[(*numc)];
X (*numc)++;
X k = Endetacol[(*numc)];
X theta = 1 / (double) Eta_value[k - 1];
X Eta_value[k - 1] = theta;
X for (i = j; i < k - 1; i++)
X Eta_value[i] *= -theta;
X} /* addetacol */
X
X
Xvoid presolve(void)
X{
X int i, j, rownr;
X double theta;
X int *num, *rownum;
X
X CALLOC(num, Rows + 2, int);
X CALLOC(rownum, Rows + 2, int);
X
X if (Verbose)
X printf("presolve\n");
X Rh[0] = 0;
X for (i = 0; i <= Sum; i++)
X Lower[i] = TRUE;
X for (i = 1; i <= Rows; i++)
X Basis[i] = TRUE;
X for (i = 1; i <= Columns; i++)
X Basis[Rows + i] = FALSE;
X for (i = 0; i <= Rows; i++)
X Bas[i] = i;
X
X /* we want to maximise the objective function, so: */
X Chsign[0] = TRUE;
X
X for (i = 1; i <= Rows; i++)
X if ((Relat[i] == GE) && (Upbo[i] == INFINITE))
X Chsign[i] = TRUE;
X else
X Chsign[i] = FALSE;
X
X /* invert values in rows with Chsign[i] == TRUE */
X for (i = 1; i <= Rows; i++)
X if (Chsign[i])
X Rh[i] = -Rh[i];
X for (i = 0; i < Nonnuls; i++)
X if (Chsign[Mat[i].rownr])
X Mat[i].value = -Mat[i].value;
X
X /* handle lower bounds by transforming to 0 */
X for (i = 1; i <= Columns; i++)
X if (Lowbo[Rows + i] > 0)
X {
X theta = Lowbo[Rows + i];
X Upbo[Rows + i] -= theta;
X for (j = Cend[i - 1]; j < Cend[i]; j++)
X Rh[Mat[j].rownr] -= theta * Mat[j].value;
X }
X
X for (i = 1; i <= Rows; i++)
X {
X num[i] = 0;
X rownum[i] = 0;
X }
X for (i = 0; i < Nonnuls; i++)
X rownum[Mat[i].rownr]++;
X Rend[0] = 0;
X for (i = 1; i <= Rows; i++)
X Rend[i] = Rend[i - 1] + rownum[i];
X for (i = 1; i <= Columns; i++)
X for (j = Cend[i - 1]; j < Cend[i]; j++)
X {
X rownr = Mat[j].rownr;
X if (rownr != 0)
X {
X num[rownr]++;
X Colno[Rend[rownr - 1] + num[rownr]] = i;
X }
X }
X free((char *)num);
X free((char *)rownum);
X} /* presolve */
X
Xvoid setpivcol(short *lower,
X int varin,
X int numeta,
X double *pcol)
X{
X int i, colnr;
X double f;
X
X if (Verbose)
X printf("setpivcol\n");
X if ((*lower))
X f = 1;
X else
X f = -1;
X for (i = 0; i <= Rows; i++)
X pcol[i] = 0;
X if (varin > Rows)
X {
X colnr = varin - Rows;
X for (i = Cend[colnr - 1]; i < Cend[colnr]; i++)
X pcol[Mat[i].rownr] = Mat[i].value * f;
X pcol[0] -= Extrad * f;
X }
X else
X if ((*lower))
X pcol[varin] = 1;
X else
X pcol[varin] = -1;
X ftran(1, numeta, pcol);
X} /* setpivcol */
X
X
Xshort colprim(int *colnr,
X int *numeta,
X short *minit,
X double *drow)
X{
X int varnr, i, j;
X double f, dpiv;
X
X if (Verbose)
X printf("colprim\n");
X dpiv = -EPSD;
X (*colnr) = 0;
X if (!(*minit))
X {
X for (i = 1; i <= Rows; i++)
X drow[i] = 0;
X drow[0] = 1;
X btran((*numeta), drow);
X for (i = 1; i <= Columns; i++)
X {
X varnr = Rows + i;
X if (!Basis[varnr])
X if (Upbo[varnr] > 0)
X {
X f = 0;
X for (j = Cend[i - 1]; j < Cend[i]; j++)
X f += drow[Mat[j].rownr] * Mat[j].value;
X drow[varnr] = f;
X }
X }
X for (i = 1; i <= Sum; i++)
X if (abs(drow[i]) < EPSD)
X drow[i] = 0;
X }
X for (i = 1; i <= Sum; i++)
X if (!Basis[i])
X if (Upbo[i] > 0)
X {
X if (Lower[i])
X f = drow[i];
X else
X f = -drow[i];
X if (f < dpiv)
X {
X dpiv = f;
X (*colnr) = i;
X }
X }
X return ((*colnr) > 0);
X} /* colprim */
X
X
Xvoid minoriteration(int colnr,
X int rownr,
X int *numeta)
X{
X int i, j, k, wk, varin, varout, elnr;
X double pivot, theta;
X
X if (Verbose)
X printf("minoriteration\n");
X varin = colnr + Rows;
X elnr = Endetacol[(*numeta)];
X wk = elnr;
X (*numeta)++;
X if (Extrad != 0)
X {
X Eta_rownr[elnr] = 0;
X Eta_value[elnr] = -Extrad;
X elnr++;
X }
X for (j = Cend[colnr - 1] ; j < Cend[colnr]; j++)
X {
X k = Mat[j].rownr;
X if (k == 0 && Extrad != 0)
X Eta_value[Endetacol[(*numeta) - 1]] += Mat[j].value;
X else
X if (k != rownr)
X {
X Eta_rownr[elnr] = k;
X Eta_value[elnr] = Mat[j].value;
X elnr++;
X }
X else
X pivot = Mat[j].value;
X }
X Eta_rownr[elnr] = rownr;
X Eta_value[elnr] = 1 / (double) pivot;
X elnr++;
X theta = Rhs[rownr] / (double) pivot;
X Rhs[rownr] = theta;
X for (i = wk; i < elnr - 1; i++)
X Rhs[Eta_rownr[i]] -= theta * Eta_value[i];
X varout = Bas[rownr];
X Bas[rownr] = varin;
X Basis[varout] = FALSE;
X Basis[varin] = TRUE;
X for (i = wk; i < elnr - 1; i++)
X Eta_value[i] /= - (double) pivot;
X /* printf("minoriteration: new etaindex: %d\n", elnr); */
X Endetacol[(*numeta)] = elnr;
X} /* minoriteration */
X
X
Xvoid rhsmincol(double *theta,
X int rownr,
X int varin,
X int numeta)
X{
X int i, j, k, varout;
X double f;
X
X if (Verbose)
X printf("rhsmincol\n");
X if (rownr > Rows + 1)
X {
X fprintf(stderr, "Error: rhsmincol called with rownr: %d, Rows: %d\n",
X rownr, Rows);
X }
X j = Endetacol[numeta];
X k = Endetacol[numeta + 1];
X for (i = j; i < k; i++)
X {
X f = Rhs[Eta_rownr[i]] - (*theta) * Eta_value[i];
X if (abs(f) < EPSB)
X Rhs[Eta_rownr[i]] = 0;
X else
X Rhs[Eta_rownr[i]] = f;
X }
X Rhs[rownr] = (*theta);
X varout = Bas[rownr];
X Bas[rownr] = varin;
X Basis[varout] = FALSE;
X Basis[varin] = TRUE;
X} /* rhsmincol */
X
Xvoid invert(int n,
X int *numeta,
X int *numinv)
X{
X int i, j, v, wk, numit, varnr, rownr, colnr, varin;
X double f, theta;
X double *pcol;
X short *frow;
X short *fcol;
X int *rownum, *col, *row;
X int *colnum;
X short *lbas;
X double *rhs1;
X
X CALLOC(rownum, Rows + 1, int);
X CALLOC(col, Rows + 1, int);
X CALLOC(row, Rows + 1, int);
X CALLOC(pcol, Rows + 2, double);
X CALLOC(frow, Rows + 2, short);
X CALLOC(fcol, Columns + 1, short);
X CALLOC(colnum, Columns + 2, int);
X CALLOC(lbas, Sum +2, short);
X CALLOC(rhs1, Rows + 2, double);
X
X if (Verbose)
X printf("a: invert:%d %d %12f\n", (*numeta), Endetacol[(*numeta)], Rhs[0]);
X
X for(i = 0; i <= Sum + 1; i++)
X lbas[i] = Basis[i];
X for(i = 0; i <= Rows + 1; i++)
X rhs1[i] = Rhs[i];
X for (i = 0; i <= Rows; i++)
X frow[i] = TRUE;
X for (i = 0; i < n; i++)
X fcol[i] = FALSE;
X for (i = 0; i < Rows; i++)
X rownum[i] = 0;
X for (i = 1; i <= n; i++)
X colnum[i] = 0;
X for (i = 0; i <= Rows; i++)
X if (Bas[i] > Rows)
X fcol[Bas[i] - Rows - 1] = TRUE;
X else
X frow[Bas[i]] = FALSE;
X for (i = 1; i <= Rows; i++)
X if (frow[i])
X for (j = Rend[i - 1] + 1; j <= Rend[i]; j++)
X {
X wk = Colno[j];
X if (fcol[wk - 1])
X {
X colnum[wk]++;
X rownum[i - 1]++;
X }
X }
X for (i = 1; i <= Rows; i++)
X Bas[i] = i;
X for (i = 1; i <= Rows; i++)
X Basis[i] = TRUE;
X for (i = 1; i <= n; i++)
X Basis[i + Rows] = FALSE;
X for (i = 0; i <= Rows; i++)
X Rhs[i] = Rh[i];
X for (i = 1; i <= n; i++)
X {
X varnr = Rows + i;
X if (!Lower[varnr])
X {
X theta = Upbo[varnr];
X for (j = Cend[i - 1]; j < Cend[i]; j++)
X Rhs[Mat[j].rownr] -= theta * Mat[j].value;
X }
X }
X for (i = 1; i <= Rows; i++)
X if (!Lower[i])
X Rhs[i] -= Upbo[i];
X (*numeta) = 0;
X v = 0;
X rownr = 0;
X (*numinv) = 0;
X numit = 0;
X while (v < Rows)
X {
X rownr++;
X if (rownr > Rows)
X rownr = 1;
X v++;
X if (rownum[rownr - 1] == 1)
X if (frow[rownr])
X {
X v = 0;
X j = Rend[rownr - 1] + 1;
X while (!(fcol[Colno[j] - 1]))
X j++;
X colnr = Colno[j];
X fcol[colnr - 1] = FALSE;
X colnum[colnr] = 0;
X for (j = Cend[colnr - 1]; j < Cend[colnr]; j++)
X if (frow[Mat[j].rownr])
X rownum[Mat[j].rownr - 1]--;
X frow[rownr] = FALSE;
X minoriteration(colnr, rownr, numeta);
X }
X }
X v = 0;
X colnr = 0;
X while (v < n)
X {
X colnr++;
X if (colnr > n)
X colnr = 1;
X v++;
X if (colnum[colnr] == 1)
X if (fcol[colnr - 1])
X {
X v = 0;
X j = Cend[colnr - 1] + 1;
X while (!(frow[Mat[j - 1].rownr]))
X j++;
X rownr = Mat[j - 1].rownr;
X frow[rownr] = FALSE;
X rownum[rownr - 1] = 0;
X for (j = Rend[rownr - 1] + 1; j <= Rend[rownr]; j++)
X if (fcol[Colno[j] - 1])
X colnum[Colno[j]]--;
X fcol[colnr - 1] = FALSE;
X numit++;
X col[numit - 1] = colnr;
X row[numit - 1] = rownr;
X }
X }
X for (j = 1; j <= n; j++)
X if (fcol[j - 1])
X {
X fcol[j - 1] = FALSE;
X setpivcol(&Lower[Rows + j], j + Rows, (*numeta), pcol);
X rownr = 1;
X while (!(frow[rownr] && pcol[rownr]))
X rownr++; /* this sometimes sets rownr to Rows + 2 and makes */
X /* rhsmincol crash. MB */
X frow[rownr] = FALSE;
X condensecol(rownr, (*numeta), pcol);
X theta = Rhs[rownr] / (double) pcol[rownr];
X rhsmincol(&theta, rownr, Rows + j, (*numeta));
X addetacol(numeta);
X }
X for (i = numit - 1; i >= 0; i--)
X {
X colnr = col[i];
X rownr = row[i];
X varin = colnr + Rows;
X for (j = 0; j <= Rows; j++)
X pcol[j] = 0;
X for (j = Cend[colnr - 1]; j < Cend[colnr]; j++)
X pcol[Mat[j].rownr] = Mat[j].value;
X pcol[0] -= Extrad;
X condensecol(rownr, (*numeta), pcol);
X theta = Rhs[rownr] / (double) pcol[rownr];
X rhsmincol(&theta, rownr, varin, (*numeta));
X addetacol(numeta);
X }
X for (i = 1; i <= Rows; i++)
X if (abs(Rhs[i]) < EPSB)
X Rhs[i] = 0;
X f = 0;
X for (i = 1; i <= Rows; i++)
X if (Rhs[i] < 0)
X f += Rhs[i];
X else
X if (Rhs[i] > Upbo[Bas[i]])
X f = f + Upbo[Bas[i]] - Rhs[i];
X if (Verbose)
X printf("b: invert:%d %d %12f %12f\n", (*numeta), Endetacol[(*numeta)],
X Rhs[0], f);
X
X free((char *)rownum);
X free((char *)col);
X free((char *)row);
X free((char *)pcol);
X free((char *)frow);
X free((char *)fcol);
X free((char *)colnum);
X free((char *)lbas);
X free((char *)rhs1);
X} /* invert */
X
X
Xshort rowprim(int *rownr,
X double *theta,
X double *pcol)
X{
X int i;
X double f, quot;
X
X if (Verbose)
X printf("rowprim\n");
X (*rownr) = 0;
X (*theta) = INFINITE;
X for (i = 1; i <= Rows; i++)
X {
X f = pcol[i];
X if (f != 0)
X {
X quot = (*theta) * 10;
X if (f > 0)
X quot = Rhs[i] / (double) f;
X else
X if (Upbo[Bas[i]] < INFINITE)
X quot = (Rhs[i] - Upbo[Bas[i]]) / (double) f;
X if (quot < (*theta))
X {
X (*theta) = quot;
X (*rownr) = i;
X }
X }
X }
X return ((*rownr) > 0);
X} /* rowprim */
X
X
X
Xvoid iteration(int *numc,
X int *rownr,
X int *varin,
X double *theta,
X double *up,
X short *minit,
X short *low,
X short *primair,
X short *smotes,
X int *iter,
X int *numinv)
X{
X int i, k, varout;
X double f, pivot;
X
X if (Verbose)
X printf("iteration\n");
X (*iter)++;
X (*minit) = (*theta) > (*up);
X if ((*minit))
X {
X (*theta) = (*up);
X (*low) = !(*low);
X }
X k = Endetacol[(*numc) + 1];
X pivot = Eta_value[k - 1];
X for (i = Endetacol[(*numc)]; i < k; i++)
X {
X f = Rhs[Eta_rownr[i]] - (*theta) * Eta_value[i];
X if (abs(f) < EPSB)
X Rhs[Eta_rownr[i]] = 0;
X else
X Rhs[Eta_rownr[i]] = f;
X }
X if (!(*minit))
X {
X Rhs[(*rownr)] = (*theta);
X varout = Bas[(*rownr)];
X Bas[(*rownr)] = (*varin);
X Basis[varout] = FALSE;
X Basis[(*varin)] = TRUE;
X if ((*primair) && pivot < 0)
X Lower[varout] = FALSE;
X if (!(*low) && (*up) != INFINITE)
X {
X (*low) = TRUE;
X Rhs[(*rownr)] = (*up) - Rhs[(*rownr)];
X for (i = Endetacol[(*numc)]; i < k; i++)
X Eta_value[i] = -Eta_value[i];
X }
X addetacol(numc);
X (*numinv)++;
X }
X if ((*smotes) && Verbose)
X {
X printf("Iteration %d: ", (*iter));
X if ((*minit))
X {
X printf("%4d", (*varin));
X if (Lower[(*varin)])
X printf(" ");
X else
X printf(" u");
X printf("%22c", ' ');
X }
X else
X {
X printf("%4d", (*varin));
X if (Lower[(*varin)])
X printf(" ");
X else
X printf(" u");
X printf("%4d", varout);
X if (Lower[varout])
X printf(" ");
X else
X printf(" u");
X printf("%12.5f", pivot);
X }
X printf(" %12f", Pcol[0]);
X if (!(*primair))
X {
X f = 0;
X for (i = 1; i <= Rows; i++)
X if (Rhs[i] < 0)
X f += Rhs[i];
X else
X if (Rhs[i] > Upbo[Bas[i]])
X f = f + Upbo[Bas[i]] - Rhs[i];
X printf(" %12f", f);
X }
X else
X printf(" %12f", Rhs[0]);
X printf("\n");
X }
X} /* iteration */
X
X
Xint solvelp(void)
X{
X int i, j, iter, varnr;
X double f, theta;
X short primair, doiter, doinvert, smotes;
X double *drow, *prow;
X int numinv, invnum;
X int numeta;
X short artif, minit;
X int colnr, rownr;
X
X CALLOC(drow, Sum + 2, double);
X CALLOC(prow, Sum + 2, double);
X
X numinv = 0;
X invnum = 50; /* number of iterations between inversions */
X numeta = 0;
X iter = 0;
X minit = FALSE;
X primair = TRUE;
X i = 0;
X smotes = TRUE;
X
X while (i != Rows && primair)
X {
X i++;
X primair = Rhs[i] >= 0 && Rhs[i] <= Upbo[Bas[i]];
X }
X
X if (!primair)
X {
X drow[0] = 1;
X for (i = 1; i <= Rows; i++)
X drow[i] = 0;
X Extrad = 0;
X for (i = 1; i <= Columns; i++)
X {
X varnr = Rows + i;
X drow[varnr] = 0;
X for (j = Cend[i - 1]; j < Cend[i]; j++)
X if (drow[Mat[j].rownr] != 0)
X drow[varnr] += drow[Mat[j].rownr] * Mat[j].value;
X if (drow[varnr] < Extrad)
X Extrad = drow[varnr];
X }
X artif = Extrad < -EPSD;
X }
X else
X {
X artif = FALSE;
X Extrad = 0;
X }
X if (Verbose)
X printf("artificial:%12.5f\n", Extrad);
X minit = FALSE;
X do {
X doiter = FALSE;
X doinvert = FALSE;
X if (primair)
X {
X rownr = 0;
X if (colprim(&colnr, &numeta, &minit, drow))
X {
X setpivcol(&Lower[colnr], colnr, numeta, Pcol);
X if (rowprim(&rownr, &theta, Pcol))
X {
X doiter = TRUE;
X condensecol(rownr, numeta, Pcol);
X }
X }
X }
X else
X {
X if (!minit)
X rowdual(&rownr);
X if (!doinvert)
X if (rownr > 0)
X {
X if (coldual(&numeta, &rownr, &colnr, &minit, prow, drow))
X {
X doiter = TRUE;
X setpivcol(&Lower[colnr], colnr, numeta, Pcol);
X condensecol(rownr, numeta, Pcol);
X f = Rhs[rownr] - Upbo[Bas[rownr]];
X if (f > 0)
X {
X theta = f / (double) Pcol[rownr];
X if (theta <= Upbo[colnr])
X Lower[Bas[rownr]] = !Lower[Bas[rownr]];
X }
X else /* getting div by zero here ... MB */
X theta = Rhs[rownr] / (double) Pcol[rownr];
X }
X }
X else
X if (artif)
X {
X primair = TRUE;
X artif = FALSE;
X printf("primal\n");
X doinvert = TRUE;
X Extrad = 0;
X }
X else
X colnr = 0;
X }
X if (doiter)
X iteration(&numeta, &rownr, &colnr, &theta, &Upbo[colnr],
X &minit, &Lower[colnr], &primair, &smotes, &iter, &numinv);
X if (numinv >= invnum)
X doinvert = TRUE;
X if (doinvert)
X invert(Columns, &numeta, &numinv);
X } while (rownr && colnr || doinvert);
X
X free((char *)drow);
X free((char *)prow);
X
X return(rownr || colnr); /* return nonzero if anything was wrong */
X} /* solvelp */
X
X
Xvoid construct_solution(double *sol)
X{
X int i, j, basi;
X double f;
X
X /* zero all results of rows */
X bzero((char *)sol, (Rows + 1) * sizeof(double));
X
X for (i = Rows + 1; i <= Sum; i++)
X sol[i] = Lowbo[i];
X for (i = 1; i <= Rows; i++)
X {
X basi = Bas[i];
X if (basi > Rows)
X sol[basi] += Rhs[i];
X }
X for (i = Rows + 1; i <= Sum; i++)
X if (!Basis[i] && !Lower[i])
X sol[i] += Upbo[i];
X for (j = 1; j <= Columns; j++)
X {
X f = sol[Rows + j];
X if (f != 0)
X for (i = Cend[j - 1]; i < Cend[j]; i++)
X sol[Mat[i].rownr] += f * Mat[i].value;
X }
X
X for (i = 0; i <= Rows; i++)
X {
X if (abs(sol[i]) < EPSB)
X sol[i] = 0;
X else
X if (Chsign[i])
X sol[i] = -sol[i];
X }
X} /* construct_solution */
X
X
Xint solve(double *upbo,
X double *lowbo)
X{
X int i, failure, notint;
X intrec *ptr;
X
X Level++;
X debug_print("starting solve");
X
X /* make fresh copies of Mat, Rh, Upbo, Lowbo, as solving changes them */
X bcopy((char*)Orig_mat, (char *)Mat, (Nonnuls + 1) * sizeof(matrec));
X bcopy((char*)Orig_rh, (char *)Rh, (Rows + 2) * sizeof(double));
X bcopy((char*)upbo, (char *)Upbo, (Sum + 2) * sizeof(double));
X bcopy((char*)lowbo, (char *)Lowbo, (Sum + 2) * sizeof(double));
X
X Rend[0] = 0;
X presolve();
X
X for (i = 0; i <= Rows; i++)
X Rhs[i] = Rh[i];
X
X Endetacol[0] = 0;
X failure = solvelp();
X
X if(failure)
X debug_print("this problem has no solution");
X
X if(!failure) /* there is a solution */
X {
X construct_solution(Solution);
X
X debug_print("a solution was found");
X debug_print_solution(Solution);
X
X /* if this solution is worse than the best sofar, this branch must die */
X if(Solution[0] <= Best_solution[0])
X {
X debug_print("but it was worse than the best sofar, discarded");
X Level--;
X return(TRUE);
X }
X
X /* check if solution contains enough ints */
X notint = 0;
X for(ptr = First_int; !notint && ptr; ptr = ptr->next)
X if(!is_int(Solution[ptr->varnr]))
X notint = ptr->varnr;
X
X if(notint) /* there is at least one value not yet int */
X {
X /* set up two new problems */
X double *new_upbo, *new_lowbo;
X double new_bound;
X int res1, res2;
X
X /* allocate room for them */
X MALLOC(new_upbo, Sum + 2, double);
X MALLOC(new_lowbo, Sum + 2, double);
X
X bcopy((char *)upbo, (char *)new_upbo, (Sum + 2) * sizeof(double));
X bcopy((char *)lowbo, (char *)new_lowbo, (Sum + 2) * sizeof(double));
X
X debug_print("not enough ints. Selecting var %s, val: %10.3g",
X Names[notint], Solution[notint]);
X
X new_bound = floor(Solution[notint]);
X new_upbo[notint] = new_bound;
X
X debug_print("starting first subproblem with bounds:");
X debug_print_bounds(new_upbo, lowbo);
X
X res1 = solve(new_upbo, lowbo);
X
X new_bound += 1;
X new_lowbo[notint] = new_bound;
X
X debug_print("starting second subproblem with bounds:");
X debug_print_bounds(upbo, new_lowbo);
X
X res2 = solve(upbo, new_lowbo);
X
X if(res1 && res2) /* both failed */
X failure = TRUE;
X else
X failure = FALSE;
X
X FREE(new_upbo);
X FREE(new_lowbo);
X }
X else /* all required values are int */
X {
X debug_print("--> Valid solution found");
X if(Solution[0] > Best_solution[0]) /* Current solution better */
X {
X debug_print("the best sofar. Prev: %10.3g, New: %10.3g",
X Best_solution[0], Solution[0]);
X bcopy((char *)Solution, (char *)Best_solution,
X (Sum + 1) * sizeof(double));
X if(Show_results)
X {
X fprintf(stderr, "Intermediate solution:\n");
X print_solution(stderr, Best_solution);
X }
X }
X }
X }
X Level--;
X return(failure);
X} /* solve */
END_OF_lp_solve/solve.c
if test 18876 -ne `wc -c <lp_solve/solve.c`; then
echo shar: \"lp_solve/solve.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
echo shar: End of archive 2 \(of 2\).
cp /dev/null ark2isdone
MISSING=""
for I in 1 2 ; do
if test ! -f ark${I}isdone ; then
MISSING="${MISSING} ${I}"
fi
done
if test "${MISSING}" = "" ; then
echo You have unpacked both archives.
rm -f ark[1-9]isdone
else
echo You still need to unpack the following archives:
echo " " ${MISSING}
fi
## End of shell archive.
exit 0
exit 0 # Just in case...
--
Andrew Patrick acting as Comp.Sources.Reviewed Moderator
Department of Communications, Ottawa, CANADA
csr@calvin.dgbt.doc.CA
Thread Navigation
This is a paginated view of messages in the thread with full content displayed inline.
Messages are displayed in chronological order, with the original post highlighted in green.
Use pagination controls to navigate through all messages in large threads.
Back to All Threads