/[ascend]/trunk/lp_solve_2.0/lpkit.c
ViewVC logotype

Annotation of /trunk/lp_solve_2.0/lpkit.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (hide annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (21 years, 1 month ago) by aw0a
File MIME type: text/x-csrc
File size: 44686 byte(s)
Setting up web subdirectory in repository
1 aw0a 1 #include "lpkit.h"
2     #include "lpglob.h"
3     #include <stdlib.h>
4     #include <stdarg.h>
5     #include <stdio.h>
6    
7     /* Globals */
8     lprec *Lp=NULL; /* pointer to active problem */
9     int Rows;
10     int Columns;
11     int Sum;
12     int Non_zeros;
13     int Level;
14     matrec *Mat;
15     int *Col_no;
16     int *Col_end;
17     int *Row_end;
18     REAL *Orig_rh;
19     REAL *Rh;
20     REAL *Rhs;
21     short *Must_be_int;
22     REAL *Orig_upbo;
23     REAL *Orig_lowbo;
24     REAL *Upbo;
25     REAL *Lowbo;
26     int *Bas;
27     short *Basis;
28     short *Lower;
29     int Eta_alloc;
30     int Eta_size;
31     REAL *Eta_value;
32     int *Eta_row_nr;
33     int *Eta_col_end;
34     int Num_inv;
35     REAL *Solution;
36     REAL *Best_solution;
37     REAL Infinite;
38     REAL Epsilon;
39     REAL Epsb;
40     REAL Epsd;
41     REAL Epsel;
42    
43     REAL TREJ;
44     REAL TINV;
45    
46     short Maximise;
47     short Floor_first;
48     REAL Extrad;
49    
50     int Warn_count; /* used in CHECK version of rounding macro */
51    
52     void error(char *format, ...)
53     {
54     va_list ap;
55     va_start(ap, format);
56     vfprintf(stderr, format, ap);
57     fputc('\n', stderr);
58     va_end(ap);
59     exit(FAIL);
60     }
61    
62     lprec *make_lp(int rows, int columns)
63     {
64     lprec *newlp;
65     int i, sum;
66    
67     sum=rows+columns;
68     if(rows < 0 || columns < 0)
69     error("rows < 0 or columns < 0");
70     CALLOC(newlp, 1, lprec);
71    
72     strcpy(newlp->lp_name, "unnamed");
73     newlp->active=FALSE;
74     newlp->verbose=FALSE;
75     newlp->print_duals=FALSE;
76     newlp->print_sol=FALSE;
77     newlp->debug=FALSE;
78     newlp->print_at_invert=FALSE;
79     newlp->trace=FALSE;
80    
81     newlp->rows=rows;
82     newlp->columns=columns;
83     newlp->sum=sum;
84     newlp->rows_alloc=rows;
85     newlp->columns_alloc=columns;
86     newlp->sum_alloc=sum;
87     newlp->names_used=FALSE;
88    
89     newlp->obj_bound=DEF_INFINITE;
90     newlp->infinite=DEF_INFINITE;
91     newlp->epsilon=DEF_EPSILON;
92     newlp->epsb=DEF_EPSB;
93     newlp->epsd=DEF_EPSD;
94     newlp->epsel=DEF_EPSEL;
95     newlp->non_zeros=0;
96     newlp->mat_alloc=1;
97     CALLOC(newlp->mat, newlp->mat_alloc, matrec);
98     CALLOC(newlp->col_no, newlp->mat_alloc, int);
99     CALLOC(newlp->col_end, columns + 1, int);
100     CALLOC(newlp->row_end, rows + 1, int);
101     newlp->row_end_valid=FALSE;
102     CALLOC(newlp->orig_rh, rows + 1, REAL);
103     CALLOC(newlp->rh, rows + 1, REAL);
104     CALLOC(newlp->rhs, rows + 1, REAL);
105     CALLOC(newlp->must_be_int, sum + 1, short);
106     for(i = 0; i <= sum; i++)
107     newlp->must_be_int[i]=FALSE;
108     CALLOC(newlp->orig_upbo, sum + 1, REAL);
109     for(i = 0; i <= sum; i++)
110     newlp->orig_upbo[i]=newlp->infinite;
111     CALLOC(newlp->upbo, sum + 1, REAL);
112     CALLOC(newlp->orig_lowbo, sum + 1, REAL);
113     CALLOC(newlp->lowbo, sum + 1, REAL);
114    
115     newlp->basis_valid=TRUE;
116     CALLOC(newlp->bas, rows+1, int);
117     CALLOC(newlp->basis, sum + 1, short);
118     CALLOC(newlp->lower, sum + 1, short);
119     for(i = 0; i <= rows; i++)
120     {
121     newlp->bas[i]=i;
122     newlp->basis[i]=TRUE;
123     }
124     for(i = rows + 1; i <= sum; i++)
125     newlp->basis[i]=FALSE;
126     for(i = 0 ; i <= sum; i++)
127     newlp->lower[i]=TRUE;
128    
129     newlp->eta_valid=TRUE;
130     newlp->eta_size=0;
131     newlp->eta_alloc=10000;
132     newlp->max_num_inv=DEFNUMINV;
133    
134     newlp->nr_lagrange=0;
135    
136     CALLOC(newlp->eta_value, newlp->eta_alloc, REAL);
137     CALLOC(newlp->eta_row_nr, newlp->eta_alloc, int);
138     CALLOC(newlp->eta_col_end, newlp->rows_alloc + newlp->max_num_inv, int);
139    
140     newlp->bb_rule=FIRST_NI;
141     newlp->break_at_int=FALSE;
142     newlp->break_value=0;
143    
144     newlp->iter=0;
145     newlp->total_iter=0;
146     CALLOC(newlp->solution, sum + 1, REAL);
147     CALLOC(newlp->best_solution, sum + 1, REAL);
148     CALLOC(newlp->duals, rows + 1, REAL);
149    
150     newlp->maximise = FALSE;
151     newlp->floor_first = TRUE;
152    
153     newlp->scaling_used = FALSE;
154     newlp->columns_scaled = FALSE;
155    
156     CALLOC(newlp->ch_sign, rows + 1, short);
157    
158     for(i = 0; i <= rows; i++)
159     newlp->ch_sign[i] = FALSE;
160    
161     newlp->valid = FALSE;
162    
163     return(newlp);
164     }
165    
166     void delete_lp(lprec *lp)
167     {
168     int i;
169    
170     if(lp->active)
171     Lp=NULL;
172     if(lp->names_used)
173     {
174     free(lp->row_name);
175     free(lp->col_name);
176     }
177     free(lp->mat);
178     free(lp->col_no);
179     free(lp->col_end);
180     free(lp->row_end);
181     free(lp->orig_rh);
182     free(lp->rh);
183     free(lp->rhs);
184     free(lp->must_be_int);
185     free(lp->orig_upbo);
186     free(lp->orig_lowbo);
187     free(lp->upbo);
188     free(lp->lowbo);
189     free(lp->bas);
190     free(lp->basis);
191     free(lp->lower);
192     free(lp->eta_value);
193     free(lp->eta_row_nr);
194     free(lp->eta_col_end);
195     free(lp->solution);
196     free(lp->best_solution);
197     free(lp->duals);
198     free(lp->ch_sign);
199     if(lp->scaling_used)
200     free(lp->scale);
201     if(lp->nr_lagrange>0)
202     {
203     free(lp->lag_rhs);
204     free(lp->lambda);
205     free(lp->lag_con_type);
206     for(i=0; i < lp->nr_lagrange; i++)
207     free(lp->lag_row[i]);
208     free(lp->lag_row);
209     }
210    
211     free(lp);
212     lp=NULL;
213     }
214    
215     lprec *copy_lp(lprec *lp)
216     {
217     lprec *newlp;
218     int i, rowsplus, colsplus, sumplus;
219    
220     rowsplus=lp->rows_alloc+1;
221     colsplus=lp->columns_alloc+1;
222     sumplus=lp->sum_alloc+1;
223    
224     MALLOCCPY(newlp, lp, 1, lprec); /* copy all non pointers */
225    
226     newlp->active=FALSE;
227    
228     if(newlp->names_used)
229     {
230     MALLOCCPY(newlp->col_name, lp->col_name, colsplus, nstring);
231     MALLOCCPY(newlp->row_name, lp->row_name, rowsplus, nstring);
232     }
233    
234     MALLOCCPY(newlp->mat, lp->mat, newlp->mat_alloc, matrec);
235     MALLOCCPY(newlp->col_end, lp->col_end, colsplus, int);
236     MALLOCCPY(newlp->col_no, lp->col_no, newlp->mat_alloc, int);
237     MALLOCCPY(newlp->row_end, lp->row_end, rowsplus, int);
238     MALLOCCPY(newlp->orig_rh, lp->orig_rh, rowsplus, REAL);
239     MALLOCCPY(newlp->rh, lp->rh, rowsplus, REAL);
240     MALLOCCPY(newlp->rhs, lp->rhs, rowsplus, REAL);
241     MALLOCCPY(newlp->must_be_int, lp->must_be_int, sumplus, short);
242     MALLOCCPY(newlp->orig_upbo, lp->orig_upbo, sumplus, REAL);
243     MALLOCCPY(newlp->orig_lowbo, lp->orig_lowbo, sumplus, REAL);
244     MALLOCCPY(newlp->upbo, lp->upbo, sumplus, REAL);
245     MALLOCCPY(newlp->lowbo, lp->lowbo, sumplus, REAL);
246     MALLOCCPY(newlp->bas, lp->bas, rowsplus, int);
247     MALLOCCPY(newlp->basis, lp->basis, sumplus, short);
248     MALLOCCPY(newlp->lower, lp->lower, sumplus, short);
249     MALLOCCPY(newlp->eta_value, lp->eta_value, lp->eta_alloc, REAL);
250     MALLOCCPY(newlp->eta_row_nr, lp->eta_row_nr, lp->eta_alloc, int);
251     MALLOCCPY(newlp->eta_col_end, lp->eta_col_end,
252     lp->rows_alloc + lp->max_num_inv, int);
253     MALLOCCPY(newlp->solution, lp->solution, sumplus, REAL);
254     MALLOCCPY(newlp->best_solution, lp->best_solution, sumplus, REAL);
255     MALLOCCPY(newlp->duals, lp->duals, rowsplus, REAL);
256     MALLOCCPY(newlp->ch_sign, lp->ch_sign, rowsplus, short);
257    
258     if(newlp->scaling_used)
259     MALLOCCPY(newlp->scale, lp->scale, sumplus, REAL);
260    
261     if(newlp->nr_lagrange > 0)
262     {
263     MALLOCCPY(newlp->lag_rhs, lp->lag_rhs, newlp->nr_lagrange, REAL);
264     MALLOCCPY(newlp->lambda, lp->lambda, newlp->nr_lagrange, REAL);
265     MALLOCCPY(newlp->lag_con_type, lp->lag_con_type, newlp->nr_lagrange,
266     short);
267     MALLOC(newlp->lag_row, newlp->nr_lagrange, REAL*);
268     for(i = 0; i < newlp->nr_lagrange; i++)
269     MALLOCCPY(newlp->lag_row[i], lp->lag_row[i], colsplus, REAL);
270     }
271     return(newlp);
272     }
273    
274     void inc_mat_space(lprec *lp, int maxextra)
275     {
276     if(lp->non_zeros + maxextra >= lp->mat_alloc)
277     {
278     lp->mat_alloc = lp->non_zeros + maxextra;
279     REALLOC(lp->mat, lp->mat_alloc, matrec);
280     REALLOC(lp->col_no, lp->mat_alloc, int);
281     if (lp->active)
282     {
283     Mat=lp->mat;
284     Col_no=lp->col_no;
285     }
286     }
287     }
288    
289     void inc_row_space(lprec *lp)
290     {
291     if(lp->rows > lp->rows_alloc)
292     {
293     lp->rows_alloc=lp->rows+10;
294     lp->sum_alloc=lp->rows_alloc+lp->columns_alloc;
295     REALLOC(lp->orig_rh, lp->rows_alloc + 1, REAL);
296     REALLOC(lp->rh, lp->rows_alloc + 1, REAL);
297     REALLOC(lp->rhs, lp->rows_alloc + 1, REAL);
298     REALLOC(lp->orig_upbo, lp->sum_alloc + 1, REAL);
299     REALLOC(lp->upbo, lp->sum_alloc + 1, REAL);
300     REALLOC(lp->orig_lowbo, lp->sum_alloc + 1, REAL);
301     REALLOC(lp->lowbo, lp->sum_alloc + 1, REAL);
302     REALLOC(lp->solution, lp->sum_alloc + 1, REAL);
303     REALLOC(lp->best_solution, lp->sum_alloc + 1, REAL);
304     REALLOC(lp->row_end, lp->rows_alloc + 1, int);
305     REALLOC(lp->basis, lp->sum_alloc + 1, short);
306     REALLOC(lp->lower, lp->sum_alloc + 1, short);
307     REALLOC(lp->must_be_int, lp->sum_alloc + 1, short);
308     REALLOC(lp->bas, lp->rows_alloc + 1, int);
309     REALLOC(lp->duals, lp->rows_alloc + 1, REAL);
310     REALLOC(lp->ch_sign, lp->rows_alloc + 1, short);
311     REALLOC(lp->eta_col_end, lp->rows_alloc + lp->max_num_inv, int);
312     if(lp->names_used)
313     REALLOC(lp->row_name, lp->rows_alloc + 1, nstring);
314     if(lp->scaling_used)
315     REALLOC(lp->scale, lp->sum_alloc + 1, REAL);
316     if(lp->active)
317     set_globals(lp);
318     }
319     }
320    
321     void inc_col_space(lprec *lp)
322     {
323     if(lp->columns >= lp->columns_alloc)
324     {
325     lp->columns_alloc=lp->columns+10;
326     lp->sum_alloc=lp->rows_alloc+lp->columns_alloc;
327     REALLOC(lp->must_be_int, lp->sum_alloc + 1, short);
328     REALLOC(lp->orig_upbo, lp->sum_alloc + 1, REAL);
329     REALLOC(lp->upbo, lp->sum_alloc + 1, REAL);
330     REALLOC(lp->orig_lowbo, lp->sum_alloc + 1, REAL);
331     REALLOC(lp->lowbo, lp->sum_alloc + 1, REAL);
332     REALLOC(lp->solution, lp->sum_alloc + 1, REAL);
333     REALLOC(lp->best_solution, lp->sum_alloc + 1, REAL);
334     REALLOC(lp->basis, lp->sum_alloc + 1, short);
335     REALLOC(lp->lower, lp->sum_alloc + 1, short);
336     if(lp->names_used)
337     REALLOC(lp->col_name, lp->columns_alloc + 1, nstring);
338     if(lp->scaling_used)
339     REALLOC(lp->scale, lp->sum_alloc + 1, REAL);
340     REALLOC(lp->col_end, lp->columns_alloc + 1, int);
341     if(lp->active)
342     set_globals(lp);
343     }
344     }
345    
346     void set_mat(lprec *lp, int Row, int Column, REAL Value)
347     {
348     int elmnr, lastelm, i;
349    
350     if(Row > lp->rows || Row < 0)
351     error("Row out of range");
352     if(Column > lp->columns || Column < 1)
353     error("Column out of range");
354     if(lp->scaling_used)
355     Value *= lp->scale[Row] * lp->scale[lp->rows + Column];
356    
357     if(TRUE /*abs(Value) > lp->epsilon*/)
358     {
359     if (lp->basis[Column] == TRUE && Row > 0)
360     lp->basis_valid = FALSE;
361     lp->eta_valid = FALSE;
362     elmnr = lp->col_end[Column-1];
363     while((elmnr < lp->col_end[Column]) ?
364     (lp->mat[elmnr].row_nr != Row) : FALSE)
365     elmnr++;
366    
367     if((elmnr != lp->col_end[Column]) ?
368     (lp->mat[elmnr].row_nr == Row) : FALSE )
369     if (lp->scaling_used)
370     {
371     if(lp->ch_sign[Row])
372     lp->mat[elmnr].value =
373     -Value * lp->scale[Row] * lp->scale[Column];
374     else
375     lp->mat[elmnr].value =
376     Value * lp->scale[Row] * lp->scale[Column];
377     }
378     else
379     {
380     if(lp->ch_sign[Row])
381     lp->mat[elmnr].value = -Value;
382     else
383     lp->mat[elmnr].value = Value;
384     }
385     else
386     {
387     /* check if more space is needed for matrix */
388     inc_mat_space(lp,1);
389    
390     /* Shift the matrix */
391     lastelm=lp->non_zeros;
392     for(i = lastelm; i > elmnr ; i--)
393     lp->mat[i]=lp->mat[i-1];
394     for(i = Column; i <= lp->columns; i++)
395     lp->col_end[i]++;
396    
397     /* Set new element */
398     lp->mat[elmnr].row_nr=Row;
399    
400     if (lp->scaling_used)
401     {
402     if(lp->ch_sign[Row])
403     lp->mat[elmnr].value=-Value*lp->scale[Row]*lp->scale[Column];
404     else
405     lp->mat[elmnr].value=Value*lp->scale[Row]*lp->scale[Column];
406     }
407     else
408     {
409     if(lp->ch_sign[Row])
410     lp->mat[elmnr].value=-Value;
411     else
412     lp->mat[elmnr].value=Value;
413     }
414    
415     lp->row_end_valid=FALSE;
416    
417     lp->non_zeros++;
418     if (lp->active)
419     Non_zeros=lp->non_zeros;
420     }
421     }
422     }
423    
424     void set_obj_fn(lprec *lp, REAL *row)
425     {
426     int i;
427     for(i = 1; i <= lp->columns; i++)
428     set_mat(lp, 0, i, row[i]);
429     }
430    
431     void str_set_obj_fn(lprec *lp, char *row)
432     {
433     int i;
434     REAL *arow;
435     char *p, *newp;
436     CALLOC(arow, lp->columns + 1, REAL);
437     p = row;
438     for(i = 1; i <= lp->columns; i++)
439     {
440     arow[i] = (REAL) strtod(p, &newp);
441     if(p==newp)
442     error("Bad string in str_set_obj_fn");
443     else
444     p=newp;
445     }
446     set_obj_fn(lp, arow);
447     free(arow);
448     }
449    
450    
451     void add_constraint(lprec *lp, REAL *row, short constr_type, REAL rh)
452     {
453     matrec *newmat;
454     int i, j;
455     int elmnr;
456     int stcol;
457     int *addtoo;
458     int *num,*rownum, row_nr;
459     REAL theta;
460    
461     MALLOC(addtoo, lp->columns + 1, int)
462    
463     for(i = 1; i <= lp->columns; i++)
464     if(row[i]!=0)
465     {
466     addtoo[i]=TRUE;
467     lp->non_zeros++;
468     }
469     else
470     addtoo[i]=FALSE;
471    
472     MALLOC(newmat, lp->non_zeros, matrec);
473     inc_mat_space(lp, 0);
474     lp->rows++;
475     lp->sum++;
476     inc_row_space(lp);
477    
478     if(lp->scaling_used)
479     {
480     /* shift scale */
481     for(i=lp->sum; i > lp->rows; i--)
482     lp->scale[i]=lp->scale[i-1];
483     lp->scale[lp->rows]=1;
484     }
485    
486     if(lp->names_used)
487     sprintf(lp->row_name[lp->rows], "r_%d", lp->rows);
488    
489     if(lp->scaling_used && lp->columns_scaled)
490     for(i = 1; i <= lp->columns; i++)
491     row[i] *= lp->scale[lp->rows+i];
492    
493     if(constr_type==GE)
494     lp->ch_sign[lp->rows] = TRUE;
495     else
496     lp->ch_sign[lp->rows] = FALSE;
497    
498     elmnr = 0;
499     stcol = 0;
500     for(i = 1; i <= lp->columns; i++)
501     {
502     for(j = stcol; j < lp->col_end[i]; j++)
503     {
504     newmat[elmnr].row_nr=lp->mat[j].row_nr;
505     newmat[elmnr].value=lp->mat[j].value;
506     elmnr++;
507     }
508     if(addtoo[i])
509     {
510     if(lp->ch_sign[lp->rows])
511     newmat[elmnr].value = -row[i];
512     else
513     newmat[elmnr].value = row[i];
514     newmat[elmnr].row_nr = lp->rows;
515     elmnr++;
516     }
517     stcol=lp->col_end[i];
518     lp->col_end[i]=elmnr;
519     }
520    
521     memcpy(lp->mat, newmat, lp->non_zeros*sizeof(matrec));
522    
523     free(newmat);
524     free(addtoo);
525    
526     for(i=lp->sum ; i > lp->rows; i--)
527     {
528     lp->orig_upbo[i]=lp->orig_upbo[i-1];
529     lp->orig_lowbo[i]=lp->orig_lowbo[i-1];
530     lp->basis[i]=lp->basis[i-1];
531     lp->lower[i]=lp->lower[i-1];
532     lp->must_be_int[i]=lp->must_be_int[i-1];
533     }
534    
535     for(i= 1 ; i <= lp->rows; i++)
536     if(lp->bas[i] >= lp->rows)
537     lp->bas[i]++;
538    
539     if(constr_type==LE || constr_type==GE)
540     {
541     lp->orig_upbo[lp->rows]=lp->infinite;
542     }
543     else if(constr_type==EQ)
544     {
545     lp->orig_upbo[lp->rows]=0;
546     }
547     else
548     {
549     fprintf(stderr, "Wrong constraint type\n");
550     exit(FAIL);
551     }
552    
553     lp->orig_lowbo[lp->rows]=0;
554    
555     if(constr_type==GE && rh != 0)
556     lp->orig_rh[lp->rows]=-rh;
557     else
558     lp->orig_rh[lp->rows]=rh;
559    
560     lp->row_end_valid=FALSE;
561    
562     lp->bas[lp->rows]=lp->rows;
563     lp->basis[lp->rows]=TRUE;
564     lp->lower[lp->rows]=TRUE;
565    
566     if(lp->active)
567     set_globals(lp);
568     lp->eta_valid=FALSE;
569     }
570    
571     void str_add_constraint(lprec *lp,
572     char *row_string,
573     short constr_type,
574     REAL rh)
575     {
576     int i;
577     REAL *aRow;
578     char *p, *newp;
579     CALLOC(aRow, lp->columns + 1, REAL);
580     p = row_string;
581    
582     for(i = 1; i <= lp->columns; i++)
583     {
584     aRow[i] = (REAL) strtod(p, &newp);
585     if(p==newp)
586     error("Bad string in str_add_constr");
587     else
588     p=newp;
589     }
590     add_constraint(lp, aRow, constr_type, rh);
591     free(aRow);
592     }
593    
594     void del_constraint(lprec *lp, int del_row)
595     {
596     int i, j;
597     unsigned elmnr;
598     int startcol;
599     int row_nr;
600    
601     if(del_row<1 || del_row>lp->rows)
602     {
603     fprintf(stderr, "There is no constraint nr. %d\n", del_row);
604     exit(FAIL);
605     }
606    
607     elmnr=0;
608     startcol=0;
609    
610     for(i = 1; i <= lp->columns; i++)
611     {
612     for(j=startcol; j < lp->col_end[i]; j++)
613     {
614     if(lp->mat[j].row_nr!=del_row)
615     {
616     lp->mat[elmnr]=lp->mat[j];
617     if(lp->mat[elmnr].row_nr > del_row)
618     lp->mat[elmnr].row_nr--;
619     elmnr++;
620     }
621     else
622     lp->non_zeros--;
623     }
624     startcol=lp->col_end[i];
625     lp->col_end[i]=elmnr;
626     }
627     for(i = del_row; i < lp->rows; i++)
628     {
629     lp->orig_rh[i] = lp->orig_rh[i+1];
630     lp->ch_sign[i] = lp->ch_sign[i+1];
631     lp->bas[i] = lp->bas[i+1];
632     if(lp->names_used)
633     strcpy(lp->row_name[i], lp->row_name[i+1]);
634     }
635     for(i = 1; i < lp->rows; i++)
636     if(lp->bas[i] > del_row)
637     lp->bas[i]--;
638    
639     for(i=del_row; i < lp->sum; i++)
640     {
641     lp->lower[i]=lp->lower[i+1];
642     lp->basis[i]=lp->basis[i+1];
643     lp->orig_upbo[i]=lp->orig_upbo[i+1];
644     lp->orig_lowbo[i]=lp->orig_lowbo[i+1];
645     lp->must_be_int[i]=lp->must_be_int[i+1];
646     if(lp->scaling_used)
647     lp->scale[i]=lp->scale[i+1];
648     }
649    
650     lp->rows--;
651     lp->sum--;
652    
653     lp->row_end_valid=FALSE;
654    
655     if(lp->active)
656     set_globals(lp);
657     lp->eta_valid=FALSE;
658     lp->basis_valid=FALSE;
659     }
660    
661     void add_lag_con(lprec *lp, REAL *row, short con_type, REAL rhs)
662     {
663     int i;
664     REAL sign;
665     if(con_type == LE || con_type == EQ)
666     sign=1;
667     else if(con_type == GE)
668     sign=-1;
669     else
670     error("con_type not implemented\n");
671    
672     lp->nr_lagrange++;
673     if(lp->nr_lagrange==1)
674     {
675     CALLOC(lp->lag_row, lp->nr_lagrange, REAL*);
676     CALLOC(lp->lag_rhs, lp->nr_lagrange, REAL);
677     CALLOC(lp->lambda, lp->nr_lagrange, REAL);
678     CALLOC(lp->lag_con_type, lp->nr_lagrange, short);
679     }
680     else
681     {
682     REALLOC(lp->lag_row, lp->nr_lagrange, REAL*);
683     REALLOC(lp->lag_rhs, lp->nr_lagrange, REAL);
684     REALLOC(lp->lambda, lp->nr_lagrange, REAL);
685     REALLOC(lp->lag_con_type, lp->nr_lagrange, short);
686     }
687     CALLOC(lp->lag_row[lp->nr_lagrange-1], lp->columns+1, REAL);
688     lp->lag_rhs[lp->nr_lagrange-1]=rhs * sign;
689     for( i=1; i <= lp->columns; i++)
690     lp->lag_row[lp->nr_lagrange-1][i]=row[i] * sign;
691     lp->lambda[lp->nr_lagrange-1]=0;
692     lp->lag_con_type[lp->nr_lagrange-1]=(con_type == EQ);
693     }
694    
695     void str_add_lag_con(lprec *lp, char *row, short con_type, REAL rhs)
696     {
697     int i;
698     REAL *a_row;
699     char *p, *new_p;
700     CALLOC(a_row, lp->columns + 1, REAL);
701     p = row;
702    
703     for(i = 1; i <= lp->columns; i++)
704     {
705     a_row[i] = (REAL) strtod(p, &new_p);
706     if(p==new_p)
707     error("Bad string in str_add_lag_con");
708     else
709     p=new_p;
710     }
711     add_lag_con(lp, a_row, con_type, rhs);
712     free(a_row);
713     }
714    
715    
716     void add_column(lprec *lp, REAL *column)
717     {
718     int i, j, elmnr, row_nr, *num, *rownum;
719    
720     lp->columns++;
721     lp->sum++;
722     inc_col_space(lp);
723     inc_mat_space(lp, lp->rows+1);
724    
725     if(lp->scaling_used)
726     {
727     for(i = 0; i <= lp->rows; i++)
728     column[i]*=lp->scale[i];
729     lp->scale[lp->sum]=1;
730     }
731    
732     elmnr=lp->col_end[lp->columns-1];
733     for(i = 0 ; i <= lp->rows ; i++)
734     if(column[i] != 0)
735     {
736     lp->mat[elmnr].row_nr=i;
737     if(lp->ch_sign[i])
738     lp->mat[elmnr].value=-column[i];
739     else
740     lp->mat[elmnr].value=column[i];
741     lp->non_zeros++;
742     elmnr++;
743     }
744     lp->col_end[lp->columns]=elmnr;
745     lp->orig_lowbo[lp->sum]=0;
746     lp->orig_upbo[lp->sum]=lp->infinite;
747     lp->lower[lp->sum]=TRUE;
748     lp->basis[lp->sum]=FALSE;
749     lp->must_be_int[lp->sum]=FALSE;
750     if(lp->names_used)
751     sprintf(lp->col_name[lp->columns], "var_%d", lp->columns);
752    
753    
754     lp->row_end_valid=FALSE;
755    
756     if(lp->active)
757     {
758     Sum=lp->sum;
759     Columns=lp->columns;
760     Non_zeros=lp->non_zeros;
761     }
762     }
763    
764     void str_add_column(lprec *lp, char *col_string)
765     {
766     int i;
767     REAL *aCol;
768     char *p, *newp;
769     CALLOC(aCol, lp->rows + 1, REAL);
770     p = col_string;
771    
772     for(i = 0; i <= lp->rows; i++)
773     {
774     aCol[i] = (REAL) strtod(p, &newp);
775     if(p==newp)
776     error("Bad string in str_add_column");
777     else
778     p=newp;
779     }
780     add_column(lp, aCol);
781     free(aCol);
782     }
783    
784     void del_column(lprec *lp, int column)
785     {
786     int i, j, from_elm, to_elm, elm_in_col;
787     if(column > lp->columns || column < 1)
788     error("Column out of range in del_column");
789     for(i = 1; i <= lp->rows; i++)
790     {
791     if(lp->bas[i]==lp->rows+column)
792     lp->basis_valid=FALSE;
793     else if(lp->bas[i] > lp->rows+column)
794     lp->bas[i]--;
795     }
796     for(i = lp->rows+column; i < lp->sum; i++)
797     {
798     if(lp->names_used)
799     strcpy(lp->col_name[i-lp->rows], lp->col_name[i-lp->rows+1]);
800     lp->must_be_int[i]=lp->must_be_int[i+1];
801     lp->orig_upbo[i]=lp->orig_upbo[i+1];
802     lp->orig_lowbo[i]=lp->orig_lowbo[i+1];
803     lp->upbo[i]=lp->upbo[i+1];
804     lp->lowbo[i]=lp->lowbo[i+1];
805     lp->basis[i]=lp->basis[i+1];
806     lp->lower[i]=lp->lower[i+1];
807     if(lp->scaling_used)
808     lp->scale[i]=lp->scale[i+1];
809     }
810     for(i = 0; i < lp->nr_lagrange; i++)
811     for(j = column; j <= lp->columns; j++)
812     lp->lag_row[i][j]=lp->lag_row[i][j+1];
813     to_elm=lp->col_end[column-1];
814     from_elm=lp->col_end[column];
815     elm_in_col=from_elm-to_elm;
816     for(i = from_elm; i < lp->non_zeros; i++)
817     {
818     lp->mat[to_elm]=lp->mat[i];
819     to_elm++;
820     }
821     for(i = column; i < lp->columns; i++)
822     lp->col_end[i]=lp->col_end[i+1]-elm_in_col;
823     lp->non_zeros -= elm_in_col;
824     lp->row_end_valid=FALSE;
825     lp->eta_valid=FALSE;
826    
827     lp->sum--;
828     lp->columns--;
829     if(lp->active)
830     set_globals(lp);
831     }
832    
833     void set_upbo(lprec *lp, int column, REAL value)
834     {
835     if(column > lp->columns || column < 1)
836     error("Column out of range");
837     if(lp->scaling_used)
838     value /= lp->scale[lp->rows + column];
839     if(value < lp->orig_lowbo[lp->rows + column])
840     error("Upperbound must be >= lowerbound");
841     lp->eta_valid = FALSE;
842     lp->orig_upbo[lp->rows+column] = value;
843     }
844    
845     void set_lowbo(lprec *lp, int column, REAL value)
846     {
847     if(column > lp->columns || column < 1)
848     error("Column out of range");
849     if(lp->scaling_used)
850     value /= lp->scale[lp->rows + column];
851     if(value > lp->orig_upbo[lp->rows + column])
852     error("Upperbound must be >= lowerbound");
853     lp->eta_valid = FALSE;
854     lp->orig_lowbo[lp->rows+column] = value;
855     }
856    
857     void set_int(lprec *lp, int column, short must_be_int)
858     {
859     if(column > lp->columns || column < 1)
860     error("Column out of range");
861     lp->must_be_int[lp->rows+column]=must_be_int;
862     if(must_be_int==TRUE)
863     if(lp->columns_scaled)
864     unscale_columns(lp);
865     }
866    
867     void set_rh(lprec *lp, int row, REAL value)
868     {
869     if(row > lp->rows || row < 0)
870     error("Row out of Range");
871     if(row == 0) /* setting of RHS of OF not meaningful */
872     {
873     fprintf(stderr,
874     "Warning: attempt to set RHS of objective function, ignored\n");
875     return;
876     }
877     if(lp->scaling_used)
878     if(lp->ch_sign[row])
879     lp->orig_rh[row] = -value * lp->scale[row];
880     else
881     lp->orig_rh[row] = value * lp->scale[row];
882     else
883     if(lp->ch_sign[row])
884     lp->orig_rh[row] = -value;
885     else
886     lp->orig_rh[row] = value;
887     lp->eta_valid = FALSE;
888     }
889    
890     void set_rh_vec(lprec *lp, REAL *rh)
891     {
892     int i;
893     if(lp->scaling_used)
894     for(i = 1; i <= lp->rows; i++)
895     if(lp->ch_sign[i])
896     lp->orig_rh[i]=-rh[i]*lp->scale[i];
897     else
898     lp->orig_rh[i]=rh[i]*lp->scale[i];
899     else
900     for(i=1; i <= lp->rows; i++)
901     if(lp->ch_sign[i])
902     lp->orig_rh[i]=-rh[i];
903     else
904     lp->orig_rh[i]=rh[i];
905     lp->eta_valid=FALSE;
906     }
907    
908     void str_set_rh_vec(lprec *lp, char *rh_string)
909     {
910     int i;
911     REAL *newrh;
912     char *p, *newp;
913     CALLOC(newrh, lp->rows + 1, REAL);
914     p = rh_string;
915    
916     for(i = 1; i <= lp->rows; i++)
917     {
918     newrh[i] = (REAL) strtod(p, &newp);
919     if(p==newp)
920     error("Bad string in str_set_rh_vec");
921     else
922     p=newp;
923     }
924     set_rh_vec(lp, newrh);
925     free(newrh);
926     }
927    
928    
929     void set_maxim(lprec *lp)
930     {
931     int i;
932     if(lp->maximise==FALSE)
933     {
934     for(i = 0; i < lp->non_zeros; i++)
935     if(lp->mat[i].row_nr==0)
936     lp->mat[i].value*=-1;
937     lp->eta_valid=FALSE;
938     }
939     lp->maximise=TRUE;
940     lp->ch_sign[0]=TRUE;
941     if(lp->active)
942     Maximise=TRUE;
943     }
944    
945     void set_minim(lprec *lp)
946     {
947     int i;
948     if(lp->maximise==TRUE)
949     {
950     for(i = 0; i < lp->non_zeros; i++)
951     if(lp->mat[i].row_nr==0)
952     lp->mat[i].value = -lp->mat[i].value;
953     lp->eta_valid=FALSE;
954     }
955     lp->maximise=FALSE;
956     lp->ch_sign[0]=FALSE;
957     if(lp->active)
958     Maximise=FALSE;
959     }
960    
961     void set_constr_type(lprec *lp, int row, short con_type)
962     {
963     int i;
964     if(row > lp->rows || row < 1)
965     error("Row out of Range");
966     if(con_type==EQ)
967     {
968     lp->orig_upbo[row]=0;
969     lp->basis_valid=FALSE;
970     if(lp->ch_sign[row])
971     {
972     for(i = 0; i < lp->non_zeros; i++)
973     if(lp->mat[i].row_nr==row)
974     lp->mat[i].value*=-1;
975     lp->eta_valid=FALSE;
976     lp->ch_sign[row]=FALSE;
977     if(lp->orig_rh[row]!=0)
978     lp->orig_rh[row]*=-1;
979     }
980     }
981     else if(con_type==LE)
982     {
983     lp->orig_upbo[row]=lp->infinite;
984     lp->basis_valid=FALSE;
985     if(lp->ch_sign[row])
986     {
987     for(i = 0; i < lp->non_zeros; i++)
988     if(lp->mat[i].row_nr==row)
989     lp->mat[i].value*=-1;
990     lp->eta_valid=FALSE;
991     lp->ch_sign[row]=FALSE;
992     if(lp->orig_rh[row]!=0)
993     lp->orig_rh[row]*=-1;
994     }
995     }
996     else if(con_type==GE)
997     {
998     lp->orig_upbo[row]=lp->infinite;
999     lp->basis_valid=FALSE;
1000     if(!lp->ch_sign[row])
1001     {
1002     for(i = 0; i < lp->non_zeros; i++)
1003     if(lp->mat[i].row_nr==row)
1004     lp->mat[i].value*=-1;
1005     lp->eta_valid=FALSE;
1006     lp->ch_sign[row]=TRUE;
1007     if(lp->orig_rh[row]!=0)
1008     lp->orig_rh[row]*=-1;
1009     }
1010     }
1011     else
1012     error("Constraint type not (yet) implemented");
1013     }
1014    
1015     REAL mat_elm(lprec *lp, int row, int column)
1016     {
1017     REAL value;
1018     int elmnr;
1019     if(row < 0 || row > lp->rows)
1020     error("Row out of range in mat_elm");
1021     if(column < 1 || column > lp->columns)
1022     error("Column out of range in mat_elm");
1023     value=0;
1024     elmnr=lp->col_end[column-1];
1025     while(lp->mat[elmnr].row_nr != row && elmnr < lp->col_end[column])
1026     elmnr++;
1027     if(elmnr != lp->col_end[column])
1028     {
1029     value = lp->mat[elmnr].value;
1030     if(lp->ch_sign[row])
1031     value = -value;
1032     if(lp->scaling_used)
1033     value /= lp->scale[row] * lp->scale[lp->rows + column];
1034     }
1035     return(value);
1036     }
1037    
1038    
1039     void get_row(lprec *lp, int row_nr, REAL *row)
1040     {
1041     int i, j;
1042     REAL RowMult;
1043     if(row_nr <0 || row_nr > lp->rows)
1044     error("Row nr. out of range in get_row");
1045     for(i = 1; i <= lp->columns; i++)
1046     {
1047     row[i]=0;
1048     for(j=lp->col_end[i-1]; j < lp->col_end[i]; j++)
1049     if(lp->mat[j].row_nr==row_nr)
1050     row[i]=lp->mat[j].value;
1051     if(lp->scaling_used)
1052     row[i] /= lp->scale[lp->rows+i] * lp->scale[row_nr];
1053     }
1054     if(lp->ch_sign[row_nr])
1055     for(i=0; i <= lp->columns; i++)
1056     if(row[i]!=0)
1057     row[i] = -row[i];
1058     }
1059    
1060     void get_column(lprec *lp, int col_nr, REAL *column)
1061     {
1062     int i;
1063    
1064     if(col_nr < 1 || col_nr > lp->columns)
1065     error("Col. nr. out of range in get_column");
1066     for(i = 0; i <= lp->rows; i++)
1067     column[i]=0;
1068     for(i = lp->col_end[col_nr-1]; i < lp->col_end[col_nr]; i++)
1069     column[lp->mat[i].row_nr]=lp->mat[i].value;
1070     for(i = 0; i <= lp->rows; i++)
1071     if(column[i] !=0)
1072     {
1073     if(lp->ch_sign[i])
1074     column[i]*=-1;
1075     if(lp->scaling_used)
1076     column[i]/=(lp->scale[i] * lp->scale[lp->rows+col_nr]);
1077     }
1078     }
1079    
1080     void get_reduced_costs(lprec *lp, REAL *rc)
1081     {
1082     int varnr, i, j;
1083     REAL f;
1084    
1085     if(!lp->basis_valid)
1086     error("Not a valid basis in get_reduced_costs");
1087     set_globals(lp);
1088     if(!lp->eta_valid)
1089     invert();
1090     for(i = 1; i <= lp->sum; i++)
1091     rc[i] = 0;
1092     rc[0] = 1;
1093     btran(rc);
1094     for(i = 1; i <= lp->columns; i++)
1095     {
1096     varnr = lp->rows + i;
1097     if(!Basis[varnr])
1098     if(Upbo[varnr] > 0)
1099     {
1100     f = 0;
1101     for(j = Col_end[i - 1]; j < Col_end[i]; j++)
1102     f += rc[Mat[j].row_nr] * Mat[j].value;
1103     rc[varnr] = f;
1104     }
1105     }
1106     for(i = 1; i <= Sum; i++)
1107     my_round(rc[i], Epsd);
1108     }
1109    
1110     short is_feasible(lprec *lp, REAL *values)
1111     {
1112     int i, elmnr;
1113     REAL *this_rhs;
1114     REAL dist;
1115    
1116     if(lp->scaling_used)
1117     {
1118     for(i = lp->rows+1; i <= lp->sum; i++)
1119     if( values[i - lp->rows] < lp->orig_lowbo[i]*lp->scale[i]
1120     || values[i - lp->rows] > lp->orig_upbo[i]*lp->scale[i])
1121     return(FALSE);
1122     }
1123     else
1124     {
1125     for(i = lp->rows+1; i <= lp->sum; i++)
1126     if( values[i - lp->rows] < lp->orig_lowbo[i]
1127     || values[i - lp->rows] > lp->orig_upbo[i])
1128     return(FALSE);
1129     }
1130     CALLOC(this_rhs, lp->rows+1, REAL)
1131     for(i = 1; i <= lp->columns; i++)
1132     for(elmnr = lp->col_end[i - 1]; elmnr < lp->col_end[i]; elmnr++)
1133     this_rhs[lp->mat[elmnr].row_nr] += lp->mat[elmnr].value * values[i];
1134     for(i = 1; i <= lp->rows; i++)
1135     {
1136     dist = lp->orig_rh[i] - this_rhs[i];
1137     my_round(dist, 0.001)
1138     if((lp->orig_upbo[i] == 0 && dist != 0) || dist < 0)
1139     {
1140     free(this_rhs);
1141     return(FALSE);
1142     }
1143     }
1144     free(this_rhs);
1145     return(TRUE);
1146     }
1147    
1148     short column_in_lp(lprec *lp, REAL *testcolumn)
1149     {
1150     int i, j;
1151     short ident;
1152     REAL value;
1153    
1154     if(lp->scaling_used)
1155     for(i = 1; i <= lp->columns; i++)
1156     {
1157     ident = TRUE;
1158     j = lp->col_end[i-1];
1159     while(ident && (j < lp->col_end[i]))
1160     {
1161     value = lp->mat[j].value;
1162     if(lp->ch_sign[lp->mat[j].row_nr])
1163     value = -value;
1164     value /= lp->scale[lp->rows+i];
1165     value /= lp->scale[lp->mat[j].row_nr];
1166     value -= testcolumn[lp->mat[j].row_nr];
1167     if(my_abs(value) > 0.001) /* should be some epsilon? MB */
1168     ident=FALSE;
1169     j++;
1170     }
1171     if(ident)
1172     return(TRUE);
1173     }
1174     else
1175     for(i = 1; i <= lp->columns; i++)
1176     {
1177     ident = TRUE;
1178     j = lp->col_end[i-1];
1179     while(ident && j < lp->col_end[i])
1180     {
1181     value = lp->mat[j].value;
1182     if(lp->ch_sign[lp->mat[j].row_nr])
1183     value *= -1;
1184     value -= testcolumn[lp->mat[j].row_nr];
1185     if( my_abs(value) > 0.001 )
1186     ident=FALSE;
1187     j++;
1188     }
1189     if(ident)
1190     return(TRUE);
1191     }
1192     return(FALSE);
1193     }
1194    
1195    
1196     void print_lp(lprec *lp)
1197     {
1198     int i, j;
1199     REAL *fatmat;
1200     CALLOC(fatmat, (lp->rows + 1) * lp->columns, REAL);
1201     for(i = 1; i <= lp->columns; i++)
1202     for(j = lp->col_end[i-1]; j < lp->col_end[i]; j++)
1203     fatmat[(i - 1) * (lp->rows + 1) + lp->mat[j].row_nr]=lp->mat[j].value;
1204    
1205     printf("problem name: %s\n", lp->lp_name);
1206     printf(" ");
1207     for(j = 1; j <= lp->columns; j++)
1208     if(lp->names_used)
1209     printf("%8s ", lp->col_name[j]);
1210     else
1211     printf("Var[%3d] ", j);
1212     if(lp->maximise)
1213     {
1214     printf("\nMaximise ");
1215     for(j = 0; j < lp->columns; j++)
1216     printf("% 8.2f ",-fatmat[j*(lp->rows+1)]);
1217     }
1218     else
1219     {
1220     printf("\nMinimize ");
1221     for(j = 0; j < lp->columns; j++)
1222     printf("% 8.2f ", fatmat[j*(lp->rows+1)]);
1223     }
1224     printf("\n");
1225     for(i = 1; i <= lp->rows; i++)
1226     {
1227     if(lp->names_used)
1228     printf("%9s ", lp->row_name[i]);
1229     else
1230     printf("Row[%3d] ", i);
1231     for(j = 0; j < lp->columns; j++)
1232     if(lp->ch_sign[i] && fatmat[j*(lp->rows+1)+i] != 0)
1233     printf("% 8.2f ",-fatmat[j*(lp->rows+1)+i]);
1234     else
1235     printf("% 8.2f ", fatmat[j*(lp->rows+1)+i]);
1236     if(lp->orig_upbo[i]==lp->infinite)
1237     if(lp->ch_sign[i])
1238     printf(">= ");
1239     else
1240     printf("<= ");
1241     else
1242     printf(" = ");
1243     if(lp->ch_sign[i])
1244     printf("% 8.2f\n",-lp->orig_rh[i]);
1245     else
1246     printf("% 8.2f\n", lp->orig_rh[i]);
1247     }
1248     printf("Type ");
1249     for(i = 1; i <= lp->columns; i++)
1250     if(lp->must_be_int[lp->rows+i]==TRUE)
1251     printf(" Int ");
1252     else
1253     printf(" Real ");
1254     printf("\nupbo ");
1255     for(i = 1; i <= lp->columns; i++)
1256     if(lp->orig_upbo[lp->rows+i]==lp->infinite)
1257     printf(" Inf ");
1258     else
1259     printf("% 8.2f ", lp->orig_upbo[lp->rows+i]);
1260     printf("\nlowbo ");
1261     for(i = 1; i <= lp->columns; i++)
1262     printf("% 8.2f ", lp->orig_lowbo[lp->rows+i]);
1263     printf("\n");
1264     for(i = 0; i < lp->nr_lagrange; i++)
1265     {
1266     printf("lag[%3d] ", i);
1267     for(j = 1; j <= lp->columns; j++)
1268     printf("% 8.2f ", lp->lag_row[i][j]);
1269     if(lp->orig_upbo[i]==lp->infinite)
1270     if(lp->lag_con_type[i] == GE)
1271     printf(">= ");
1272     else if(lp->lag_con_type[i] == LE)
1273     printf("<= ");
1274     else if(lp->lag_con_type[i] == EQ)
1275     printf(" = ");
1276     printf("% 8.2f\n", lp->lag_rhs[i]);
1277     }
1278    
1279     free(fatmat);
1280     }
1281    
1282     void set_row_name(lprec *lp, int row, nstring new_name)
1283     {
1284     int i;
1285    
1286     if(!lp->names_used)
1287     {
1288     CALLOC(lp->row_name, lp->rows_alloc + 1, nstring);
1289     CALLOC(lp->col_name, lp->columns_alloc + 1, nstring);
1290     lp->names_used=TRUE;
1291     for(i = 0; i <= lp->rows; i++)
1292     sprintf(lp->row_name[i], "r_%d", i);
1293     for(i = 1; i <= lp->columns; i++)
1294     sprintf(lp->col_name[i], "var_%d", i);
1295     }
1296     strcpy(lp->row_name[row], new_name);
1297     }
1298    
1299     void set_col_name(lprec *lp, int column, nstring new_name)
1300     {
1301     int i;
1302    
1303     if(!lp->names_used)
1304     {
1305     CALLOC(lp->row_name, lp->rows_alloc + 1, nstring);
1306     CALLOC(lp->col_name, lp->columns_alloc + 1, nstring);
1307     lp->names_used=TRUE;
1308     for(i = 0; i <= lp->rows; i++)
1309     sprintf(lp->row_name[i], "r_%d", i);
1310     for(i = 1; i <= lp->columns; i++)
1311     sprintf(lp->col_name[i], "var_%d", i);
1312     }
1313     strcpy(lp->col_name[column], new_name);
1314     }
1315    
1316     static REAL minmax_to_scale(REAL min, REAL max)
1317     {
1318     REAL scale;
1319    
1320     /* should do something sensible when min or max is 0, MB */
1321     if((min == 0) || (max == 0))
1322     return((REAL)1);
1323    
1324     scale = 1 / pow(10, (log10(min) + log10(max)) / 2);
1325     return(scale);
1326     }
1327    
1328     void unscale_columns(lprec *lp)
1329     {
1330     int i, j;
1331    
1332     /* unscale mat */
1333     for(j = 1; j <= lp->columns; j++)
1334     for(i = lp->col_end[j - 1]; i < lp->col_end[j]; i++)
1335     lp->mat[i].value /= lp->scale[lp->rows + j];
1336    
1337     /* unscale bounds as well */
1338     for(i = lp->rows + 1; i < lp->sum; i++)
1339     {
1340     if(lp->orig_lowbo[i] != 0)
1341     lp->orig_lowbo[i] *= lp->scale[i];
1342     if(lp->orig_upbo[i] != lp->infinite)
1343     lp->orig_upbo[i] *= lp->scale[i];
1344     }
1345    
1346     for(i=lp->rows+1; i<= lp->sum; i++)
1347     lp->scale[i]=1;
1348     lp->columns_scaled=FALSE;
1349     lp->eta_valid=FALSE;
1350     }
1351    
1352     void unscale(lprec *lp)
1353     {
1354     int i, j;
1355    
1356     if(lp->scaling_used)
1357     {
1358    
1359     /* unscale mat */
1360     for(j = 1; j <= lp->columns; j++)
1361     for(i = lp->col_end[j - 1]; i < lp->col_end[j]; i++)
1362     lp->mat[i].value /= lp->scale[lp->rows + j];
1363    
1364     /* unscale bounds */
1365     for(i = lp->rows + 1; i < lp->sum; i++)
1366     {
1367     if(lp->orig_lowbo[i] != 0)
1368     lp->orig_lowbo[i] *= lp->scale[i];
1369     if(lp->orig_upbo[i] != lp->infinite)
1370     lp->orig_upbo[i] *= lp->scale[i];
1371     }
1372    
1373     /* unscale the matrix */
1374     for(j = 1; j <= lp->columns; j++)
1375     for(i = lp->col_end[j-1]; i < lp->col_end[j]; i++)
1376     lp->mat[i].value /= lp->scale[lp->mat[i].row_nr];
1377    
1378     /* unscale the rhs! */
1379     for(i = 0; i <= lp->rows; i++)
1380     lp->orig_rh[i] /= lp->scale[i];
1381    
1382     free(lp->scale);
1383     lp->scaling_used=FALSE;
1384     lp->eta_valid=FALSE;
1385     }
1386     }
1387    
1388    
1389     void auto_scale(lprec *lp)
1390     {
1391     int i, j, row_nr, IntUsed;
1392     REAL *row_max, *row_min, *scalechange, absval;
1393     REAL col_max, col_min, col_scale;
1394    
1395     if(!lp->scaling_used)
1396     {
1397     MALLOC(lp->scale, lp->sum_alloc + 1, REAL);
1398     for(i=0; i <= lp->sum; i++)
1399     lp->scale[i]=1;
1400     }
1401    
1402     MALLOC(row_max, lp->rows + 1, REAL);
1403     MALLOC(row_min, lp->rows + 1, REAL);
1404     MALLOC(scalechange, lp->sum + 1, REAL);
1405    
1406     /* initialise min and max values */
1407     for(i = 0; i <= lp->rows; i++)
1408     {
1409     row_max[i]=0;
1410     row_min[i]=lp->infinite;
1411     }
1412    
1413     /* calculate min and max absolute values of rows */
1414     for(j = 1; j <= lp->columns; j++)
1415     for(i = lp->col_end[j - 1]; i < lp->col_end[j]; i++)
1416     {
1417     row_nr = lp->mat[i].row_nr;
1418     absval = my_abs(lp->mat[i].value);
1419     if(absval!=0)
1420     {
1421     row_max[row_nr] = my_max(row_max[row_nr], absval);
1422     row_min[row_nr] = my_min(row_min[row_nr], absval);
1423     }
1424     }
1425     /* calculate scale factors for rows */
1426     for(i = 0; i <= lp->rows; i++)
1427     {
1428     scalechange[i] = minmax_to_scale(row_min[i], row_max[i]);
1429     lp->scale[i] *= scalechange[i];
1430     }
1431    
1432     /* now actually scale the matrix */
1433     for(j = 1; j <= lp->columns; j++)
1434     for(i = lp->col_end[j - 1]; i < lp->col_end[j]; i++)
1435     lp->mat[i].value *= scalechange[lp->mat[i].row_nr];
1436    
1437     /* and scale the rhs! */
1438     for(i = 0; i <= lp->rows; i++)
1439     lp->orig_rh[i] *= scalechange[i];
1440    
1441     free(row_max);
1442     free(row_min);
1443    
1444     IntUsed=FALSE;
1445     i=lp->rows+1;
1446     while(!IntUsed && i <= lp->sum)
1447     {
1448     IntUsed=lp->must_be_int[i];
1449     i++;
1450     }
1451     if(!IntUsed)
1452     {
1453     REAL col_max, col_min, col_scale;
1454    
1455     /* calculate scales */
1456     for(j = 1; j <= lp->columns; j++)
1457     {
1458     col_max = 0;
1459     col_min = lp->infinite;
1460     for(i = lp->col_end[j - 1]; i < lp->col_end[j]; i++)
1461     {
1462     if(lp->mat[i].value!=0)
1463     {
1464     col_max = my_max(col_max, my_abs(lp->mat[i].value));
1465     col_min = my_min(col_min, my_abs(lp->mat[i].value));
1466     }
1467     }
1468     scalechange[lp->rows + j] = minmax_to_scale(col_min, col_max);
1469     lp->scale[lp->rows + j] *= scalechange[lp->rows + j];
1470     }
1471    
1472     /* scale mat */
1473     for(j = 1; j <= lp->columns; j++)
1474     for(i = lp->col_end[j - 1]; i < lp->col_end[j]; i++)
1475     lp->mat[i].value *= scalechange[lp->rows + j];
1476    
1477     /* scale bounds as well */
1478     for(i = lp->rows + 1; i < lp->sum; i++)
1479     {
1480     if(lp->orig_lowbo[i] != 0)
1481     lp->orig_lowbo[i] /= scalechange[i];
1482     if(lp->orig_upbo[i] != lp->infinite)
1483     lp->orig_upbo[i] /= scalechange[i];
1484     }
1485     lp->columns_scaled=TRUE;
1486     }
1487     free(scalechange);
1488     lp->scaling_used=TRUE;
1489     lp->eta_valid=FALSE;
1490     }
1491    
1492     void reset_basis(lprec *lp)
1493     {
1494     lp->basis_valid=FALSE;
1495     }
1496    
1497     void print_solution(lprec *lp)
1498     {
1499     int i;
1500    
1501     fprintf(stdout, "Value of objective function: %16.10g\n",
1502     (double)lp->best_solution[0]);
1503    
1504     /* print normal variables */
1505     for(i = 1; i <= lp->columns; i++)
1506     if(lp->names_used)
1507     fprintf(stdout, "%-10s%16.5g\n", lp->col_name[i],
1508     (double)lp->best_solution[lp->rows+i]);
1509     else
1510     fprintf(stdout, "Var [%4d] %16.5g\n", i,
1511     (double)lp->best_solution[lp->rows+i]);
1512    
1513     /* print achieved constraint values */
1514     if(lp->verbose)
1515     {
1516     fprintf(stdout, "\nActual values of the constraints:\n");
1517     for(i = 1; i <= lp->rows; i++)
1518     if(lp->names_used)
1519     fprintf(stdout, "%-10s%16.5g\n", lp->row_name[i],
1520     (double)lp->best_solution[i]);
1521     else
1522     fprintf(stdout, "%Row [%4d] %16.5g\n", i,
1523     (double)lp->best_solution[i]);
1524     }
1525    
1526     if((lp->verbose || lp->print_duals))
1527     {
1528     if(lp->max_level != 1)
1529     fprintf(stdout,
1530     "These are the duals from the node that gave the optimal solution.\n");
1531     else
1532     fprintf(stdout, "\nDual values:\n");
1533     for(i = 1; i <= lp->rows; i++)
1534     if(lp->names_used)
1535     fprintf(stdout, "%-10s%16.5g\n", lp->row_name[i],
1536     (double)lp->duals[i]);
1537     else
1538     fprintf(stdout, "Row [%4d] %16.5g\n", i, (double)lp->duals[i]);
1539     }
1540     } /* Printsolution */
1541    
1542     void write_LP(lprec *lp, FILE *output)
1543     {
1544     int i, j;
1545     REAL *row;
1546     nstring var_name;
1547    
1548     MALLOC(row, lp->columns+1, REAL);
1549     if(lp->maximise)
1550     fprintf(output, "max:");
1551     else
1552     fprintf(output, "min:");
1553    
1554     get_row(lp, 0, row);
1555     for(i = 1; i <= lp->columns; i++)
1556     if(row[i] != 0)
1557     {
1558     if(row[i] == -1)
1559     fprintf(output, " -");
1560     else if(row[i] == 1)
1561     fprintf(output, " +");
1562     else
1563     fprintf(output, " %+g ", row[i]);
1564     if(lp->names_used)
1565     fprintf(output, "%s", lp->col_name[i]);
1566     else
1567     fprintf(output, "x%d", i);
1568     }
1569     fprintf(output, ";\n");
1570    
1571     for(j = 1; j <= lp->rows; j++)
1572     {
1573     if(lp->names_used)
1574     fprintf(output, "%s:", lp->row_name[j]);
1575     get_row(lp, j, row);
1576     for(i = 1; i <= lp->columns; i++)
1577     if(row[i] != 0)
1578     {
1579     if(row[i] == -1)
1580     fprintf(output, " -");
1581     else if(row[i] == 1)
1582     fprintf(output, " +");
1583     else
1584     fprintf(output, " %+g ", row[i]);
1585     if(lp->names_used)
1586     fprintf(output, "%s", lp->col_name[i]);
1587     else
1588     fprintf(output, "x%d", i);
1589     }
1590     if(lp->orig_upbo[j] == 0)
1591     fprintf(output, " =");
1592     else if(lp->ch_sign[j])
1593     fprintf(output, " >");
1594     else
1595     fprintf(output, " <");
1596     if(lp->ch_sign[j])
1597     fprintf(output, " %g;\n",-lp->orig_rh[j]);
1598     else
1599     fprintf(output, " %g;\n", lp->orig_rh[j]);
1600     }
1601     for(i = lp->rows+1; i <= lp->sum; i++)
1602     {
1603     if(lp->orig_lowbo[i] != 0)
1604     {
1605     if(lp->names_used)
1606     fprintf(output, "%s > %g;\n", lp->col_name[i - lp->rows],
1607     lp->orig_lowbo[i]);
1608     else
1609     fprintf(output, "x%d > %g;\n", i - lp->rows,
1610     lp->orig_lowbo[i]);
1611     }
1612     if(lp->orig_upbo[i] != lp->infinite)
1613     {
1614     if(lp->names_used)
1615     fprintf(output, "%s < %g;\n", lp->col_name[i - lp->rows],
1616     lp->orig_upbo[i]);
1617     else
1618     fprintf(output, "x%d < %g;\n", i - lp->rows, lp->orig_upbo[i]);
1619     }
1620     }
1621    
1622    
1623     i=1;
1624     while(!lp->must_be_int[lp->rows+i] && i <= lp->columns)
1625     i++;
1626     if(i <= lp->columns)
1627     {
1628     if(lp->names_used)
1629     fprintf(output, "\nint %s", lp->col_name[i]);
1630     else
1631     fprintf(output, "\nint x%d", i);
1632     i++;
1633     for(; i <= lp->columns; i++)
1634     if(lp->must_be_int[lp->rows+i])
1635     if(lp->names_used)
1636     fprintf(output, ",%s", lp->col_name[i]);
1637     else
1638     fprintf(output, ", x%d", i);
1639     fprintf(output, ";\n");
1640     }
1641     free(row);
1642     }
1643    
1644    
1645    
1646    
1647     void write_MPS(lprec *lp, FILE *output)
1648     {
1649     int i, j, marker;
1650     REAL *column;
1651    
1652    
1653     MALLOC(column, lp->rows+1, REAL);
1654     marker=0;
1655     fprintf(output, "NAME %s\n", lp->lp_name);
1656     fprintf(output, "ROWS\n");
1657     for(i = 0; i <= lp->rows; i++)
1658     {
1659     if(i==0)
1660     fprintf(output, " N ");
1661     else
1662     if(lp->orig_upbo[i]==lp->infinite)
1663     if(lp->ch_sign[i])
1664     fprintf(output, " G ");
1665     else
1666     fprintf(output, " L ");
1667     else
1668     fprintf(output, " E ");
1669     if(lp->names_used)
1670     fprintf(output, "%s\n", lp->row_name[i]);
1671     else
1672     fprintf(output, "r_%d\n", i);
1673     }
1674    
1675     fprintf(output, "COLUMNS\n");
1676     j = 0;
1677     for(i = 1; i <= lp->columns; i++)
1678     {
1679     if((lp->must_be_int[i+lp->rows]) && (marker % 2)==0)
1680     {
1681     fprintf(output,
1682     " MARK%04d 'MARKER' 'INTORG'\n",
1683     marker);
1684     marker++;
1685     }
1686     if((!lp->must_be_int[i+lp->rows]) && (marker % 2)==1)
1687     {
1688     fprintf(output,
1689     " MARK%04d 'MARKER' 'INTEND'\n",
1690     marker);
1691     marker++;
1692     }
1693     get_column(lp, i, column);
1694     j=0;
1695     if(lp->maximise)
1696     {
1697     if(column[j] != 0)
1698     {
1699     if(lp->names_used)
1700     fprintf(output, " %-8s %-8s %g\n", lp->col_name[i],
1701     lp->row_name[j], -column[j]);
1702     else
1703     fprintf(output, " var_%-4d r_%-6d %g\n", i, j,
1704     -column[j]);
1705     }
1706     }
1707     else
1708     {
1709     if(column[j] != 0)
1710     {
1711     if(lp->names_used)
1712     fprintf(output, " %-8s %-8s %g\n", lp->col_name[i],
1713     lp->row_name[j], column[j]);
1714     else
1715     fprintf(output, " var_%-4d r_%-6d %g\n", i, j,
1716     column[j]);
1717     }
1718     }
1719     for(j=1; j <= lp->rows; j++)
1720     if(column[j] != 0)
1721     {
1722     if(lp->names_used)
1723     fprintf(output, " %-8s %-8s %g\n", lp->col_name[i],
1724     lp->row_name[j], column[j]);
1725     else
1726     fprintf(output, " var_%-4d r_%-6d %g\n", i, j,
1727     column[j]);
1728     }
1729     }
1730     if((marker % 2) ==1)
1731     {
1732     fprintf(output, " MARK%04d 'MARKER' 'INTEND'\n",
1733     marker);
1734     marker++;
1735     }
1736    
1737     fprintf(output, "RHS\n");
1738     for(i = 1; i <= lp->rows; i++)
1739     {
1740     if(lp->ch_sign[i])
1741     {
1742     if(lp->names_used)
1743     fprintf(output, " RHS %-8s %g\n", lp->row_name[i],
1744     (double)-lp->orig_rh[i]);
1745     else
1746     fprintf(output, " RHS r_%-6d %g\n", i,
1747     (double)-lp->orig_rh[i]);
1748     }
1749     else
1750     {
1751     if(lp->names_used)
1752     fprintf(output, " RHS %-8s %g\n", lp->row_name[i],
1753     (double)lp->orig_rh[i]);
1754     else
1755     fprintf(output, " RHS r_%-6d %g\n", i,
1756     (double)lp->orig_rh[i]);
1757     }
1758     }
1759    
1760     fprintf(output, "BOUNDS\n");
1761     if(lp->names_used)
1762     for(i = lp->rows + 1; i <= lp->sum; i++)
1763     {
1764     if(lp->orig_upbo[i] < lp->infinite)
1765     fprintf(output, " UP BND %-8s %g\n",
1766     lp->col_name[i- lp->rows], (double)lp->orig_upbo[i]);
1767     if(lp->orig_lowbo[i] != 0)
1768     fprintf(output, " LO BND %-8s %g\n",
1769     lp->col_name[i- lp->rows], (double)lp->orig_lowbo[i]);
1770     }
1771     else
1772     for(i = lp->rows + 1; i <= lp->sum; i++)
1773     {
1774     if(lp->orig_upbo[i] < lp->infinite)
1775     fprintf(output, " UP BND var_%-4d %g\n",
1776     i - lp->rows, (double)lp->orig_upbo[i]);
1777     if(lp->orig_lowbo[i] != 0)
1778     fprintf(output, " LO BND var_%-4d %g\n", i - lp->rows,
1779     (double)lp->orig_lowbo[i]);
1780     }
1781     fprintf(output, "ENDATA\n");
1782     free(column);
1783     }
1784    
1785     void print_duals(lprec *lp)
1786     {
1787     int i;
1788     for(i = 1; i <= lp->rows; i++)
1789     if(lp->names_used)
1790     fprintf(stdout, "%10s [%3d] % 10.4f\n", lp->row_name[i], i,
1791     lp->duals[i]);
1792     else
1793     fprintf(stdout, "Dual [%3d] % 10.4f\n", i, lp->duals[i]);
1794     }
1795    
1796     void print_scales(lprec *lp)
1797     {
1798     int i;
1799     if(lp->scaling_used)
1800     {
1801     for(i = 0; i <= lp->rows; i++)
1802     fprintf(stdout, "Row[%3d] scaled at % 10.6f\n", i, lp->scale[i]);
1803     for(i = 1; i <= lp->columns; i++)
1804     fprintf(stdout, "Column[%3d] scaled at % 10.6f\n", i,
1805     lp->scale[lp->rows + i]);
1806     }
1807     }

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22