1 /*
2 ** lpslink.c: function to link lpsolve with Excel, S-Plus or R.
3 */
4 /*
5 ** In addition to standard "include" files we need lpkit.h, supplied by
6 ** lpsolve. This gives definitions for (for example) the "lprec" structure.
7 */
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <math.h>
11 #include <string.h>
12 #include "lp_lib.h"
13 #include <stdio.h>
14
15 /*
16 ** In R "integers" get passed as longs, whereas in S-Plus they're ints.
17 ** So the user should turn on this flag to compile for R.
18 */
19
20 #ifdef BUILDING_FOR_R
21 #define LONG_OR_INT int
22 #else
23 #define LONG_OR_INT long
24 #endif
25
26 /*
27 ** The declaration of the link function.
28 */
29
30 void lpslink (LONG_OR_INT *direction, /* 1 for max, 0 for min */
31 LONG_OR_INT *x_count, /* Number of x's */
32 double *objective, /* Objective function */
33 LONG_OR_INT *const_count, /* Number of constraints */
34 double *constraints, /* Has extra element on front */
35 LONG_OR_INT *int_count, /* Number of integer variables */
36 LONG_OR_INT *int_vec, /* Indices of int. variables */
37 LONG_OR_INT *bin_count, /* Number of binary variables */
38 LONG_OR_INT *bin_vec, /* Indices of bin. variables */
39 LONG_OR_INT *num_bin_solns, /* # of all-bin solns to find */
40 double *obj_val, /* Objective function value */
41 double *solution, /* Result of call */
42 LONG_OR_INT *presolve, /* Value of presolve */
43 LONG_OR_INT *compute_sens, /* Want sensitivity? */
44 double *sens_coef_from, /* Sens. coef. lower limit */
45 double *sens_coef_to, /* Sens. coef. upper limit */
46 double *duals, /* Dual values */
47 double *duals_from, /* Lower limit dual values */
48 double *duals_to, /* Lower limit dual values */
49 LONG_OR_INT *scale, /* lpsolve scaling parameter */
50 LONG_OR_INT *use_dense, /* Use dense constraint matrix?*/
51 LONG_OR_INT *dense_col, /* Dense constraint columns */
52 double *dense_val, /* Dense constraint values */
53 LONG_OR_INT *dense_mat_nrow, /* Dense constraint #entries */
54 double *dense_ctr, /* Dense constraint info */
55 LONG_OR_INT *use_rw_file, /* See notes below */
56 char **rw_file, /* See notes below */
57 LONG_OR_INT *status); /* Holds return value */
58
59 /*
60 ** Some globals for calling from VBScript. The caller will call lpslink,
61 */
62 static long vb_direction;
63 static long vb_x_count;
64 static double *vb_objective;
65 static long vb_const_count;
66 static double *vb_constraints;
67 static long vb_int_count;
68 static long *vb_int_vec;
69 static double vb_obj_val;
70 static double *vb_solution;
71
72 /*
73 **************************** lps_vb_setup *****************************
74 **
75 **
76 ** lps_vb_setup: set up an lp problem (that is, allocate necessary space)
77 */
78
lps_vb_setup(long direction,long x_count,long const_count,long int_count)79 long lps_vb_setup (long direction, /* Optimization direction (1 = max) */
80 long x_count, /* Number of variables */
81 long const_count, /* Number of constraints */
82 long int_count) /* Number of integer variables */
83 {
84 long i; /* iteration variable */
85
86 /*
87 ** Set globals (which start "vb") from parameters (which do not).
88 */
89
90 vb_direction = direction;
91 vb_x_count = x_count;
92 vb_const_count = const_count;
93 vb_int_count = int_count;
94
95 /*
96 ** Allocate objective (which holds the coefficients in the objective function).
97 ** We need an extra element at the front. If malloc() fails, get out.
98 */
99 vb_objective = (double *) malloc (1 + sizeof (double) * vb_x_count);
100 if (vb_objective == (double *) NULL)
101 return (-1);
102
103 vb_objective[0] = 0.0;
104
105 /*
106 ** Allocate constraints. This, too, gets an extra entry. If the allocation
107 ** fails, exit gracefully.
108 */
109
110 vb_constraints = (double *) malloc (sizeof (double) *
111 (1 + vb_const_count * (vb_x_count + 2)));
112
113 if (vb_constraints == (double *) NULL)
114 {
115 free (vb_objective);
116 return (-1);
117 }
118
119 vb_constraints[0] = 0.0;
120
121 /*
122 ** If any variables are constrained to be integers, allocate one long
123 ** for each such variable. Quit if the allocation fails. If it succeeds,
124 ** put zeros in there. We will insert the proper numbers later.
125 */
126
127 if (vb_int_count > 0) {
128 vb_int_vec = (long *) malloc (1 + sizeof (long) * vb_int_count);
129 if (vb_int_vec == (long *) NULL)
130 {
131 free (vb_objective);
132 free (vb_constraints);
133 return (-1);
134 }
135
136 for (i = 0L; i <= vb_int_count; i++) /* there's one extra */
137 vb_int_vec[i] = 0L;
138 }
139
140 /*
141 ** Allocate space for the solution. This will hold the optimal values for each
142 ** coefficient. If the allocation fails, quit.
143 */
144
145 vb_solution = (double *) malloc (sizeof (double) * vb_x_count);
146 if (vb_solution == (double *) NULL)
147 {
148 free (vb_objective);
149 free (vb_constraints);
150 if (vb_int_count > 0)
151 free (vb_int_vec);
152 return (-1);
153 }
154 /*
155 ** Our work here is done.
156 */
157 return (0);
158 } /* end lps_vb_setup */
159
160 /***************************** lps_vb_set_element ********************/
161
lps_vb_set_element(long type,long row,long column,double value)162 long lps_vb_set_element (long type, /* Place to do the setting */
163 long row, /* Row in which to set */
164 long column, /* Column in which to set */
165 double value) /* Value to set */
166 {
167 /*
168 ** This function allows us to set an element of the lp. If "type" = 1,
169 ** we ignore column and set the "row-th" element of vb_objective to "value."
170 ** If type is 2, set the row, column-th element of constraints (allowing for
171 ** the funny layout) to value; if type = 3, set the rowth variable to integer.
172 */
173
174 if (type == 1)
175 vb_objective[row] = value;
176 if (type == 2) {
177 vb_constraints[(row - 1) * (vb_x_count + 2) + column] = value;
178 }
179 if (type == 3 && vb_int_count > 0)
180 vb_int_vec[row] = floor (value + 0.5);
181 return (1);
182 }
183
184 /***************************** lps_vb_get_element ********************/
lps_vb_get_element(long type,long row,long column)185 double lps_vb_get_element (long type, /* Place to get from */
186 long row, /* Row to get from */
187 long column) /* Column to get from (unused) */
188 {
189 /*
190 ** Get an element. If type is 1, get the objective value; if type is 2,
191 ** get the rowth element of the solution. "Column" is currently unused.
192 */
193 if (type == 1)
194 return (vb_obj_val);
195 if (type == 2)
196 return (vb_solution[row]);
197 return (0.0);
198 }
199
200 /*
201 ********************************* lps_vb_cleanup *************************
202 **
203 ** lps_vb_cleanup: free all the things allocated in lps_vb_setup.
204 */
lps_vb_cleanup(long unused)205 long lps_vb_cleanup (long unused)
206 {
207 free (vb_objective);
208 free (vb_constraints);
209 free (vb_int_vec);
210 free (vb_solution);
211 return (0);
212 }
213
214 /******************************** lpslink ************************************/
215
lpslink(LONG_OR_INT * direction,LONG_OR_INT * x_count,double * objective,LONG_OR_INT * const_count,double * constraints,LONG_OR_INT * int_count,LONG_OR_INT * int_vec,LONG_OR_INT * bin_count,LONG_OR_INT * bin_vec,LONG_OR_INT * num_bin_solns,double * obj_val,double * solution,LONG_OR_INT * presolve,LONG_OR_INT * compute_sens,double * sens_coef_from,double * sens_coef_to,double * duals,double * duals_from,double * duals_to,LONG_OR_INT * scale,LONG_OR_INT * use_dense,LONG_OR_INT * dense_col,double * dense_val,LONG_OR_INT * dense_mat_nrow,double * dense_ctr,LONG_OR_INT * use_rw_file,char ** rw_file,LONG_OR_INT * status)216 void lpslink (LONG_OR_INT *direction, /* 1 for max, 0 for min */
217 LONG_OR_INT *x_count, /* Number of x's */
218 double *objective, /* Objective function */
219 LONG_OR_INT *const_count, /* Number of constraints */
220 double *constraints, /* Has extra element on front */
221 LONG_OR_INT *int_count, /* Number of integer variables */
222 LONG_OR_INT *int_vec, /* Indices of int. variables */
223 LONG_OR_INT *bin_count, /* Number of binary variables */
224 LONG_OR_INT *bin_vec, /* Indices of bin. variables */
225 LONG_OR_INT *num_bin_solns, /* # of all-bin solns to find */
226 double *obj_val, /* Objective function value */
227 double *solution, /* Result of call */
228 LONG_OR_INT *presolve, /* Value of presolve */
229 LONG_OR_INT *compute_sens, /* Want sensitivity? */
230 double *sens_coef_from, /* Sens. coef. lower limit */
231 double *sens_coef_to, /* Sens. coef. upper limit */
232 double *duals, /* Dual values */
233 double *duals_from, /* Lower limit dual values */
234 double *duals_to, /* Lower limit dual values */
235 LONG_OR_INT *scale, /* lpsolve scaling parameter */
236 LONG_OR_INT *use_dense, /* Use dense const. mat? */
237 LONG_OR_INT *dense_col, /* Dense constraint column */
238 double *dense_val, /* Dense constraint value */
239 LONG_OR_INT *dense_mat_nrow, /* Dense constraint #entries */
240 double *dense_ctr, /* Dense constraint info */
241 LONG_OR_INT *use_rw_file, /* See notes below */
242 char **rw_file, /* See notes below */
243 LONG_OR_INT *status) /* Holds return value */
244 {
245 /*
246 ** This is the function called from the outside.
247 */
248 /*
249 ** "rw file" notes: to get around some sort of a bug in lpSolve, we allow
250 ** callers to set "use_rw_file" to TRUE and pass a file name in rw_file.
251 ** This only makes sense in the case where all the decision variables are
252 ** binary and num_bin_solns > 1. Then instead of just adding constraints to
253 ** lp and re-solving, we write the lp out to the rw_file, delete the lp, and
254 ** read the file back in every time. It's costly, but what can you do?
255 */
256
257 int i = 0, j, /* Iteration variables */
258 result; /* Holds result of calls */
259
260 int dctr_ctr, /* Holds our spot in the dense_ctr matrix */
261 dmat_ctr, /* Holds the current row of the dense_mat matrix */
262 d_num; /* Number of non-zero entries in this dense constraint */
263
264 double *const_ptr; /* Points to a constraint */
265
266 double *last_soln; /* Points to last solution (for multiple soln case) */
267 double *new_ptr; /* Point to a bit of "solution" 4 building constraint */
268 int soln_ctr; /* Which solution are we on? */
269
270 double new_rhs; /* RHS value for new constraint. */
271 LONG_OR_INT
272 new_status; /* Status for calls to "solve" after the first. */
273
274 lprec *lp; /* Structure to hold the lp */
275
276 FILE *filex_file; /* Points to rw_file once that's been opened */
277
278 /*
279 ** Make an empty lp with x_count variables. If it fails, return.
280 */
281 lp = make_lp ((int) 0, *x_count);
282
283 if (lp == (lprec *) NULL)
284 return;
285
286 set_verbose (lp, 1); /* CRITICAL */
287
288 /* Set the direction. The default is minimize, but set it anyway. */
289 if (*direction == 1)
290 set_maxim (lp);
291 else
292 set_minim (lp);
293 /*
294 ** "Objective" is a vector. Set the objective function. Return on fail.
295 */
296 result = set_obj_fn (lp, objective);
297 if (result == 0)
298 return;
299
300 set_add_rowmode (lp, TRUE); /* Put problem into "row mode" */
301
302 /*
303 ** If there are any constraints, see if they're dense or regular.
304 */
305 if ((int) *const_count > 0) {
306 if (*use_dense) {
307 /*
308 ** For dense constraints, dense_ctr holds sets of three entries (# of non-
309 ** zero entries, direction, and rhs). Once we know the number of non-zero
310 ** entries, we get them out of dense_const and use that information to fill
311 ** up a row which we dispatch with add_constraintex.
312 */
313 dctr_ctr = 0;
314 dmat_ctr = 0;
315 for (i = 0; i < (int) *const_count; i++)
316 {
317 d_num = (int) dense_ctr[dctr_ctr];
318 /*
319 ** The args. to add_constraintex are the lp, the number of non-zero entries,
320 ** a pointer to the values, an int pointer to the column numbers associated with
321 ** the values, the constraint type (<, = >) and the constraint's right side.
322 */
323 add_constraintex (lp, d_num,
324 (double *) &(dense_val[dmat_ctr]),
325 (int *) &(dense_col[dmat_ctr]),
326 (int) dense_ctr[dctr_ctr + 1], dense_ctr[dctr_ctr + 2]);
327 dctr_ctr += 3;
328 dmat_ctr += d_num;
329 }
330 /* Maybe replace this with set_row, set_rh_vec, set_constr_type? */
331 }
332 else {
333 /*
334 ** If we're using regular constaints, point "constr_ptr" at the first one.
335 */
336 const_ptr = constraints;
337 /*
338 ** Add constraints, one at a time; then move constr_ptr up.
339 */
340
341 for (i = 0; i < (int) *const_count; i++)
342 {
343 add_constraint (lp, const_ptr,
344 (short) const_ptr[(int) (*x_count) + 1],
345 const_ptr[(int) (*x_count) + 2]);
346 const_ptr += (int) *x_count + 2;
347 }
348 } /* end "else" (i.e. if these are regular, non-dense constraints.) */
349 } /* end "if there are any constraints. */
350
351 set_add_rowmode (lp, FALSE); /* Take problem out of "row mode" */
352
353 /*
354 ** Set integer and binary variables.
355 */
356 if (*int_count > 0) {
357 for (i = 0; i < (int) *int_count; i++)
358 set_int (lp, (int) (int_vec[i]), TRUE);
359 }
360 if (*bin_count > 0) {
361 for (i = 0; i < (int) *bin_count; i++)
362 set_binary (lp, (int) (bin_vec[i]), TRUE);
363 }
364
365 /*
366 ** Presolve if needed (that is, if we are going to want sensitivity in
367 ** an integer program) then solve the lp. If "status" is non-zero, return.
368 */
369
370 if (*compute_sens > 0) {
371 if (*int_count > 0)
372 set_presolve (lp, PRESOLVE_SENSDUALS, get_presolveloops (lp));
373 }
374
375
376 set_scaling (lp, *scale);
377 *status = (LONG_OR_INT) solve (lp);
378
379 if ((int) *status != 0) {
380 delete_lp (lp);
381 return;
382 }
383
384 /* Now get the sensitivities, if requested. */
385 if (*compute_sens > 0) {
386 get_sensitivity_obj (lp, sens_coef_from, sens_coef_to);
387 get_sensitivity_rhs (lp, duals, duals_from, duals_to);
388 }
389
390 /*
391 ** We've succeeded. Extract the objective function's value and
392 ** the values of the variables.
393 */
394
395 *obj_val = get_objective (lp);
396
397 get_variables (lp, solution);
398
399 /*
400 ** If this is an all-binary problem, and more than one solution has
401 ** been asked for, let's get those other solutions. We do that in
402 ** this way. First, add a constraint requiring that the objective function
403 ** not move below *obj_val, if this is a maximization problem, or above
404 ** it, if this is minimization. We need only do this once.
405 */
406 soln_ctr = 1;
407 if (*num_bin_solns > 1) {
408 add_constraint (lp, objective,
409 (*direction == 1) ? ROWTYPE_GE : ROWTYPE_LE, *obj_val);
410
411
412 /* ==================================================
413 ** ------------------ WHILE LOOP --------------------
414 ** ==================================================
415 ** Loop until new status isn't 0 or we've produced all solutions asked for.
416 */
417 new_status = 0;
418 while (new_status == 0 && soln_ctr < *num_bin_solns) {
419 /*
420 ** Point to the most recent solution and space in "solution" that we'll
421 ** use to build the new constraint.
422 */
423 last_soln = &(solution[(soln_ctr - 1) * *x_count]);
424 new_ptr = last_soln + *x_count; /* Move up *x_count elements */
425 /*
426 ** Now we need to add one new constraint. We go through every element of
427 ** the most recent solution. For every element that's a 1 in the solution,
428 ** put in a 1 in the constraint. For every 0, put a -1 in the constraint.
429 ** The rhs is the # of 1's minus 1, and the sign is <. So if there were
430 ** five variables and the last solution was (1, 1, 1, 0, 0), the new
431 ** constraint would say x1 + x2 + x3 -x4 -x5 < 2.
432 */
433 new_rhs = 0;
434 new_ptr[0] = 0;
435 for (j = 0; j < *x_count; j++) {
436 new_ptr[j + 1] = 2 * last_soln[j] - 1; /* new_ptr is one-based */
437 new_rhs += last_soln[j];
438 }
439
440 if (*use_rw_file) {
441 filex_file = fopen (*rw_file, "w");
442 write_LP (lp, filex_file);
443 delete_lp (lp);
444 fclose (filex_file);
445 filex_file = fopen (*rw_file, "r");
446 lp = read_lp (filex_file, CRITICAL, (char *) NULL);
447 fclose (filex_file);
448 }
449
450 result = add_constraint (lp, new_ptr, ROWTYPE_LE, new_rhs - 1);
451
452 set_scaling (lp, *scale);
453 new_status = (LONG_OR_INT) solve (lp);
454
455 if (new_status != 0) {
456 *num_bin_solns = soln_ctr;
457 return;
458 }
459 soln_ctr ++;
460 /*
461 ** Fetch solution into new_ptr, then transfer just *x_count entries.
462 */
463 result = get_variables (lp, new_ptr);
464
465 } /* end "while" */
466
467 *num_bin_solns = soln_ctr;
468
469 } /* end "multiple solutions */
470
471
472 /*
473 **
474 */
475
476 /*
477 ** Free up the memory and return.
478 */
479
480 delete_lp (lp);
481
482 return;
483
484 } /* end "lpslink" */
485
486 /*
487 ****************************** lp_transbig ************************
488 **
489 ** This function handles "big" transportation problem. It takes in
490 ** the number of rows and columns, the cost matrix, and the signs and
491 ** right-hand sides of the constraints and builds everything else on the
492 ** fly.
493 */
lp_transbig(LONG_OR_INT * direction,LONG_OR_INT * r_count,LONG_OR_INT * c_count,double * costs,LONG_OR_INT * r_signs,double * r_rhs,LONG_OR_INT * c_signs,double * c_rhs,double * obj_val,LONG_OR_INT * int_count,LONG_OR_INT * integers,double * solution,LONG_OR_INT * presolve,LONG_OR_INT * compute_sens,double * sens_coef_from,double * sens_coef_to,double * duals,double * duals_from,double * duals_to,LONG_OR_INT * status)494 void lp_transbig (LONG_OR_INT *direction, /* 1 for max, 0 for min */
495 LONG_OR_INT *r_count, /* Number of rows */
496 LONG_OR_INT *c_count, /* Number of columns */
497 double *costs, /* Objective function */
498 LONG_OR_INT *r_signs, /* Signs of row constraints */
499 double *r_rhs, /* RHS of row constraints */
500 LONG_OR_INT *c_signs, /* Signs of col constraints */
501 double *c_rhs, /* RHS of col constraints */
502 double *obj_val, /* Objective function value */
503 LONG_OR_INT *int_count, /* How many vars are integers?*/
504 LONG_OR_INT *integers, /* Which vars. are integer? */
505 double *solution, /* Result of call */
506 LONG_OR_INT *presolve, /* Value of presolve */
507 LONG_OR_INT *compute_sens, /* Want sensitivity? */
508 double *sens_coef_from, /* Sens. coef. lower limit */
509 double *sens_coef_to, /* Sens. coef. upper limit */
510 double *duals, /* Dual values */
511 double *duals_from, /* Lower limit dual values */
512 double *duals_to, /* Lower limit dual values */
513 LONG_OR_INT *status) /* Holds return value */
514 {
515 long i; /* Iteration variable */
516 long result; /* Holds result of calls */
517 long this_element; /* Which are we looking at? */
518 lprec *lp; /* Structure to hold the lp */
519 double *row_vals; /* Holds the values for row-type constraints */
520 int *col_inds; /* Holds locations for col-type constraints */
521 double *col_vals; /* Holds the values for col-type constraints */
522 int *row_inds; /* Holds locations for row-type constraints */
523
524 long col_ind_ctr, row_ind_ctr;
525 long rc = *r_count, cc = *c_count;
526
527 /*
528 ** Make an empty lp with r_count x c_count variables. If it fails, return.
529 */
530 lp = make_lp ((int) 0, *r_count * *c_count);
531
532 if (lp == (lprec *) NULL)
533 return;
534
535 set_verbose (lp, 1); /* CRITICAL */
536
537 set_add_rowmode (lp, TRUE);
538 /*
539 ** "Costs" is already a vector. Set the objective function. Return on fail.
540 */
541 result = set_obj_fn (lp, costs);
542 if (result == 0)
543 return;
544
545 /* Set the direction. The default is minimize, but set it anyway. */
546 if (*direction == 1)
547 set_maxim (lp);
548 else
549 set_minim (lp);
550
551 /*
552 ** Add constraints. There are r_count row-type constraints, plus c_count
553 ** col_type constraints.
554 */
555 row_vals = calloc (cc, sizeof (double));
556 col_inds = calloc (cc, sizeof (int));
557
558 for (row_ind_ctr = 0L; row_ind_ctr < rc; row_ind_ctr++)
559 {
560 for (col_ind_ctr = 0; col_ind_ctr < cc; col_ind_ctr++) {
561 row_vals[col_ind_ctr] = 1.0;
562 this_element = 1 + (col_ind_ctr * rc) + row_ind_ctr;
563 col_inds[col_ind_ctr] = this_element;
564 }
565 add_constraintex (lp, cc, row_vals, col_inds, r_signs[row_ind_ctr], r_rhs[row_ind_ctr]);
566 }
567
568 free (row_vals);
569 free (col_inds);
570
571 col_vals = calloc (rc, sizeof (double));
572 row_inds = calloc (rc, sizeof (int));
573
574 for (col_ind_ctr = 0L; col_ind_ctr < cc; col_ind_ctr++)
575 {
576 for (row_ind_ctr = 0; row_ind_ctr < rc; row_ind_ctr++) {
577 col_vals[row_ind_ctr] = 1.0;
578 this_element = 1 + row_ind_ctr + col_ind_ctr * rc;
579 row_inds[row_ind_ctr] = this_element;
580 }
581 add_constraintex (lp, rc, col_vals, row_inds, c_signs[col_ind_ctr], c_rhs[col_ind_ctr]);
582 }
583 free (col_vals);
584 free (row_inds);
585
586 set_add_rowmode (lp, FALSE);
587
588 /*
589 ** Set integers.
590 */
591 if (*int_count > 0)
592 for (i = 0; i < *int_count; i++)
593 set_int (lp, integers[i], 1); /* Variable in ith element of integers */
594
595 if (*compute_sens > 0) {
596 set_presolve (lp, PRESOLVE_SENSDUALS, 10);
597 }
598
599 *status = (LONG_OR_INT) solve (lp);
600
601 if ((int) *status != 0) {
602 return;
603 }
604
605 /* Now get the sensitivities, if requested. */
606 if (*compute_sens > 0) {
607 get_sensitivity_obj (lp, sens_coef_from, sens_coef_to);
608 get_sensitivity_rhs (lp, duals, duals_from, duals_to);
609 }
610
611 /*
612 ** We've succeeded. Extract the objective function's value and
613 ** the values of the variables.
614 */
615
616 *obj_val = get_objective (lp);
617
618 get_variables (lp, solution);
619
620 /*
621 **
622 */
623
624 /*
625 ** Free up the memory and return.
626 */
627
628 delete_lp (lp);
629 }
630