1
2 #include <string.h>
3 #include "commonlib.h"
4 #include "lp_lib.h"
5 #include "lp_report.h"
6 #include "lp_pricePSE.h"
7 #include "lp_price.h"
8
9 #if libBLAS > 0
10 #include "myblas.h"
11 #endif
12
13 #ifdef FORTIFY
14 # include "lp_fortify.h"
15 #endif
16
17 /* Simplex pricing utility module - w/interface for lp_solve v5.0+
18 -------------------------------------------------------------------------
19 Author: Kjell Eikland
20 Contact: kjell.eikland@broadpark.no
21 License terms: LGPL.
22
23 Requires: lp_lib.h, commonlib.h
24
25 Release notes:
26 v1.0.0 1 July 2004 Routines extracted from lp_lib.
27 v1.0.1 10 July 2004 Added comparison operators for determination
28 of entering and leaving variables.
29 Added routines for multiple and partial
30 pricing and made corresponding changes to
31 colprim and rowdual.
32 v1.0.2 20 August 2004 Implemented relative pivot size control in
33 rowprim and rowdual.
34 v1.1.0 15 October 2004 Added dual long step logic.
35 v1.1.1 22 October 2004 Added bound sort order to variable selections.
36 v1.2.0 24 March 2005 Completed multiple pricing logic.
37 ------------------------------------------------------------------------- */
38
39 INLINE REAL normalizeEdge(lprec *lp, int item, REAL edge, MYBOOL isdual);
40
41 /* Comparison operators for entering and leaving variables for both the primal and
42 dual simplexes. The functions compare a candidate variable with an incumbent. */
compareImprovementVar(const pricerec * current,const pricerec * candidate)43 int CMP_CALLMODEL compareImprovementVar(const pricerec *current, const pricerec *candidate)
44 {
45 register int result = COMP_PREFERNONE;
46 register lprec *lp = current->lp;
47 register REAL testvalue, margin = PREC_IMPROVEGAP;
48 int currentcolno, currentvarno = current->varno,
49 candidatecolno, candidatevarno = candidate->varno;
50 MYBOOL isdual = candidate->isdual;
51
52 if(isdual) {
53 candidatevarno = lp->var_basic[candidatevarno];
54 currentvarno = lp->var_basic[currentvarno];
55 }
56 candidatecolno = candidatevarno - lp->rows;
57 currentcolno = currentvarno - lp->rows;
58
59 /* Do pivot-based selection unless Bland's (first index) rule is active */
60 if(lp->_piv_rule_ != PRICER_FIRSTINDEX) {
61
62 MYBOOL candbetter;
63
64 /* Find the largest value - normalize in case of the dual, since
65 constraint violation is expressed as a negative number. */
66 /* Use absolute test for "small numbers", relative otherwise */
67 testvalue = candidate->pivot;
68 if(fabs(testvalue) < LIMIT_ABS_REL)
69 testvalue -= current->pivot;
70 else
71 testvalue = my_reldiff(testvalue, current->pivot);
72 if(isdual)
73 testvalue = -testvalue;
74
75 candbetter = (MYBOOL) (testvalue > 0);
76 if(candbetter) {
77 if(testvalue > margin)
78 result = COMP_PREFERCANDIDATE;
79 }
80 #if 0 /* Give more opportunity to optimize on non-primary criteria */
81 else if (testvalue < -margin)
82 #else /* Give reduced opportunity to optimize on non-primary criteria */
83 else if (testvalue < -lp->epsvalue)
84 #endif
85 result = COMP_PREFERINCUMBENT;
86
87 #ifdef UseSortOnBound
88 /* Extra selection criterion based on the variable's range;
89 variable with - DUAL: small bound out; PRIMAL: large bound in */
90 if(result == COMP_PREFERNONE) {
91 testvalue = lp->upbo[candidatevarno] - lp->upbo[currentvarno];
92 if(testvalue < -margin)
93 result = COMP_PREFERINCUMBENT;
94 else if(testvalue > margin)
95 result = COMP_PREFERCANDIDATE;
96 result = my_chsign(isdual, result);
97 }
98 #endif
99
100 #ifdef UseSortOnColumnLength
101 /* Prevent long columns from entering the basis */
102 if(result == COMP_PREFERNONE) {
103 if(candidatecolno > 0)
104 testvalue = mat_collength(lp->matA, candidatecolno) +
105 (is_obj_in_basis(lp) && (lp->obj[candidatecolno] != 0) ? 1 : 0);
106 else
107 testvalue = 1;
108 if(currentcolno > 0)
109 testvalue -= mat_collength(lp->matA, currentcolno) +
110 (is_obj_in_basis(lp) && (lp->obj[currentcolno] != 0) ? 1 : 0);
111 else
112 testvalue -= 1;
113 if(testvalue > 0)
114 result = COMP_PREFERINCUMBENT;
115 else if(testvalue < 0)
116 result = COMP_PREFERCANDIDATE;
117 result = my_chsign(isdual, result);
118 }
119 #endif
120
121 /* Select absolute best if the non-primary criteria failed to separate */
122 if((result == COMP_PREFERNONE) && candbetter) {
123 result = COMP_PREFERCANDIDATE;
124 goto Finish;
125 }
126 }
127
128 /* Final tie-breakers */
129 if(result == COMP_PREFERNONE) {
130
131 /* Add randomization tie-braker */
132 if(lp->piv_strategy & PRICE_RANDOMIZE) {
133 result = my_sign(PRICER_RANDFACT - rand_uniform(lp, 1.0));
134 if(candidatevarno < currentvarno)
135 result = -result;
136 }
137
138 /* Resolve ties via index ordinal */
139 if(result == COMP_PREFERNONE) {
140 if(candidatevarno < currentvarno)
141 result = COMP_PREFERCANDIDATE;
142 else /* if(candidatevarno > currentvarno) */
143 result = COMP_PREFERINCUMBENT;
144 if(lp->_piv_left_)
145 result = -result;
146 }
147 }
148
149 Finish:
150 return( result );
151
152 }
153
compareSubstitutionVar(const pricerec * current,const pricerec * candidate)154 int CMP_CALLMODEL compareSubstitutionVar(const pricerec *current, const pricerec *candidate)
155 {
156 register int result = COMP_PREFERNONE;
157 register lprec *lp = current->lp;
158 register REAL testvalue = candidate->theta,
159 margin = current->theta;
160 MYBOOL isdual = candidate->isdual, candbetter;
161 int currentcolno, currentvarno = current->varno,
162 candidatecolno, candidatevarno = candidate->varno;
163
164 if(!isdual) {
165 candidatevarno = lp->var_basic[candidatevarno];
166 currentvarno = lp->var_basic[currentvarno];
167 }
168 candidatecolno = candidatevarno - lp->rows;
169 currentcolno = currentvarno - lp->rows;
170
171 /* Compute the ranking test metric. */
172 if(isdual) {
173 testvalue = fabs(testvalue);
174 margin = fabs(margin);
175 }
176
177 /* Use absolute test for "small numbers", relative otherwise */
178 if(fabs(testvalue) < LIMIT_ABS_REL)
179 testvalue -= margin;
180 else
181 testvalue = my_reldiff(testvalue, margin);
182
183 /* Find if the new Theta is smaller or near equal (i.e. testvalue <= eps)
184 compared to the previous best; ties will be broken by pivot size or index
185 NB! The margin below is essential in maintaining primal/dual feasibility
186 during the primal/dual simplex, respectively. Sometimes a small
187 value prevents the selection of a suitable pivot, thereby weakening
188 the numerical stability of some models */
189 margin = PREC_SUBSTFEASGAP;
190 candbetter = (MYBOOL) (testvalue < 0);
191 if(candbetter) {
192 if(testvalue < -margin)
193 result = COMP_PREFERCANDIDATE;
194 }
195 else if(testvalue > margin)
196 result = COMP_PREFERINCUMBENT;
197
198 /* Resolve a tie */
199 if(result == COMP_PREFERNONE) {
200 REAL currentpivot = fabs(current->pivot),
201 candidatepivot = fabs(candidate->pivot);
202
203 /* Handle first index / Bland's rule specially */
204 if(lp->_piv_rule_ == PRICER_FIRSTINDEX) {
205 #if 1
206 /* Special secondary selection by pivot size (limited stability protection) */
207 margin = candidate->epspivot;
208 if((candidatepivot >= margin) && (currentpivot < margin))
209 result = COMP_PREFERCANDIDATE;
210 #endif
211 }
212
213 else {
214
215 /* General secondary selection based on pivot size */
216 #if 0
217 if(candidatepivot > MIN_STABLEPIVOT)
218 testvalue = my_reldiff(testvalue, currentpivot);
219 else
220 #endif
221 testvalue = candidatepivot - currentpivot;
222 if(testvalue > margin)
223 result = COMP_PREFERCANDIDATE;
224 else if(testvalue < -margin)
225 result = COMP_PREFERINCUMBENT;
226
227 #ifdef UseSortOnBound
228 /* Extra selection criterion based on the variable's range;
229 variable with - PRIMAL: small bound out; DUAL: large bound in */
230 if(result == COMP_PREFERNONE) {
231 testvalue = lp->upbo[candidatevarno] - lp->upbo[currentvarno];
232 if(testvalue < -margin)
233 result = COMP_PREFERCANDIDATE;
234 else if(testvalue > margin)
235 result = COMP_PREFERINCUMBENT;
236 result = my_chsign(isdual, result);
237 }
238 #endif
239
240 #ifdef UseSortOnColumnLength
241 /* Prevent long columns from entering the basis */
242 if(result == COMP_PREFERNONE) {
243 if(candidatecolno > 0)
244 testvalue = mat_collength(lp->matA, candidatecolno) +
245 (is_obj_in_basis(lp) && (lp->obj[candidatecolno] != 0) ? 1 : 0);
246 else
247 testvalue = 1;
248 if(currentcolno > 0)
249 testvalue -= mat_collength(lp->matA, currentcolno) +
250 (is_obj_in_basis(lp) && (lp->obj[currentcolno] != 0) ? 1 : 0);
251 else
252 testvalue -= 1;
253 if(testvalue > 0)
254 result = COMP_PREFERCANDIDATE;
255 else if(testvalue < 0)
256 result = COMP_PREFERINCUMBENT;
257 result = my_chsign(isdual, result);
258 }
259 #endif
260
261 }
262 }
263
264 /* Select absolute best if the non-primary criteria failed to separate */
265 if((result == COMP_PREFERNONE) && candbetter) {
266 result = COMP_PREFERCANDIDATE;
267 goto Finish;
268 }
269
270 /* Final tie-breakers */
271 if(result == COMP_PREFERNONE) {
272
273 /* Add randomization tie-braker */
274 if(lp->piv_strategy & PRICE_RANDOMIZE) {
275 result = my_sign(PRICER_RANDFACT - rand_uniform(lp, 1.0));
276 if(candidatevarno < currentvarno)
277 result = -result;
278 }
279
280 /* Resolve ties via index ordinal (also prefers slacks over user variables) */
281 if(result == COMP_PREFERNONE) {
282 if(candidatevarno < currentvarno)
283 result = COMP_PREFERCANDIDATE;
284 else /* if(candidatevarno > currentvarno) */
285 result = COMP_PREFERINCUMBENT;
286 if(lp->_piv_left_)
287 result = -result;
288 }
289 }
290
291 Finish:
292 return( result );
293 }
compareBoundFlipVar(const pricerec * current,const pricerec * candidate)294 int CMP_CALLMODEL compareBoundFlipVar(const pricerec *current, const pricerec *candidate)
295 {
296 register REAL testvalue, margin;
297 register int result = COMP_PREFERNONE;
298 register lprec *lp = current->lp;
299 MYBOOL candbetter;
300 int currentvarno = current->varno,
301 candidatevarno = candidate->varno;
302
303 if(!current->isdual) {
304 candidatevarno = lp->var_basic[candidatevarno];
305 currentvarno = lp->var_basic[currentvarno];
306 }
307
308 /* Compute the ranking test metric. */
309 testvalue = candidate->theta;
310 margin = current->theta;
311 if(candidate->isdual) {
312 testvalue = fabs(testvalue);
313 margin = fabs(margin);
314 }
315 if(fabs(margin) < LIMIT_ABS_REL)
316 testvalue -= margin;
317 else
318 testvalue = my_reldiff(testvalue, margin);
319
320 /* Find if the new Theta is smaller or near equal (i.e. testvalue <= eps)
321 compared to the previous best; ties will be broken by pivot size or index */
322 margin = PREC_SUBSTFEASGAP;
323 candbetter = (MYBOOL) (testvalue < 0);
324 if(candbetter) {
325 if(testvalue < -margin)
326 result = COMP_PREFERCANDIDATE;
327 }
328 else if(testvalue > margin)
329 result = COMP_PREFERINCUMBENT;
330
331 /* Resolve a tie */
332 if(result == COMP_PREFERNONE) {
333
334 /* Tertiary selection based on priority for large pivot sizes */
335 if(result == COMP_PREFERNONE) {
336 REAL currentpivot = fabs(current->pivot),
337 candidatepivot = fabs(candidate->pivot);
338 if(candidatepivot > currentpivot+margin)
339 result = COMP_PREFERCANDIDATE;
340 else if(candidatepivot < currentpivot-margin)
341 result = COMP_PREFERINCUMBENT;
342 }
343
344 /* Secondary selection based on priority for narrow-bounded variables */
345 if(result == COMP_PREFERNONE)
346 result = compareREAL(&(lp->upbo[currentvarno]),
347 &(lp->upbo[candidatevarno]));
348
349 }
350
351 /* Select absolute best if the non-primary criteria failed to separate */
352 if((result == COMP_PREFERNONE) && candbetter) {
353 result = COMP_PREFERCANDIDATE;
354 goto Finish;
355 }
356
357 /* Quaternary selection by index value */
358 if(result == COMP_PREFERNONE) {
359 if(candidatevarno < currentvarno)
360 result = COMP_PREFERCANDIDATE;
361 else
362 result = COMP_PREFERINCUMBENT;
363 if(lp->_piv_left_)
364 result = -result;
365 }
366
367 Finish:
368 return( result );
369 }
370
371 /* Validity operators for entering and leaving columns for both the primal and dual
372 simplex. All candidates must satisfy these tests to qualify to be allowed to be
373 a subject for the comparison functions/operators. */
validImprovementVar(pricerec * candidate)374 STATIC MYBOOL validImprovementVar(pricerec *candidate)
375 {
376 register REAL candidatepivot = fabs(candidate->pivot);
377
378 #ifdef Paranoia
379 return( (MYBOOL) ((candidate->varno > 0) && (candidatepivot > candidate->lp->epsvalue)) );
380 #else
381 return( (MYBOOL) (candidatepivot > candidate->lp->epsvalue) );
382 #endif
383 }
384
validSubstitutionVar(pricerec * candidate)385 STATIC MYBOOL validSubstitutionVar(pricerec *candidate)
386 {
387 register lprec *lp = candidate->lp;
388 register REAL theta = (candidate->isdual ? fabs(candidate->theta) : candidate->theta);
389
390 #ifdef Paranoia
391 if(candidate->varno <= 0)
392 return( FALSE );
393 else
394 #endif
395 if(fabs(candidate->pivot) >= lp->infinite)
396 return( (MYBOOL) (theta < lp->infinite) );
397 else
398 return( (MYBOOL) ((theta < lp->infinite) &&
399 (fabs(candidate->pivot) >= candidate->epspivot)) );
400 }
401
compareImprovementQS(const UNIONTYPE QSORTrec * current,const UNIONTYPE QSORTrec * candidate)402 int CMP_CALLMODEL compareImprovementQS(const UNIONTYPE QSORTrec *current, const UNIONTYPE QSORTrec *candidate)
403 {
404 return( compareImprovementVar((pricerec *) current->pvoidint2.ptr, (pricerec *) candidate->pvoidint2.ptr) );
405 }
compareSubstitutionQS(const UNIONTYPE QSORTrec * current,const UNIONTYPE QSORTrec * candidate)406 int CMP_CALLMODEL compareSubstitutionQS(const UNIONTYPE QSORTrec *current, const UNIONTYPE QSORTrec *candidate)
407 {
408 return( compareBoundFlipVar((pricerec *) current->pvoidint2.ptr, (pricerec *) candidate->pvoidint2.ptr) );
409 /* return( compareSubstitutionVar((pricerec *) current->self, (pricerec *) candidate->self) ); */
410 }
411
412 /* Function to add a valid pivot candidate into the specified list */
addCandidateVar(pricerec * candidate,multirec * multi,findCompare_func findCompare,MYBOOL allowSortedExpand)413 STATIC int addCandidateVar(pricerec *candidate, multirec *multi, findCompare_func findCompare, MYBOOL allowSortedExpand)
414 {
415 int insertpos, delta = 1;
416 pricerec *targetrec;
417
418 /* Find the insertion point (if any) */
419 if((multi->freeList[0] == 0) ||
420 (multi->sorted && allowSortedExpand) ||
421 (candidate->isdual && (multi->used == 1) && ((multi->step_last >= multi->epszero) ||
422 multi_truncatingvar(multi, ((pricerec *) (multi->sortedList[0].pvoidreal.ptr))->varno)))
423 ) {
424 UNIONTYPE QSORTrec searchTarget;
425
426 /* Make sure that the list is sorted before the search for an insertion point */
427 if((multi->freeList[0] == 0) && !multi->sorted) {
428 multi->sorted = QS_execute(multi->sortedList, multi->used, findCompare, &insertpos);
429 multi->dirty = (MYBOOL) (insertpos > 0);
430 }
431
432 /* Perform the search */
433 searchTarget.pvoidint2.ptr = (void *) candidate;
434 insertpos = sizeof(searchTarget);
435 insertpos = findIndexEx(&searchTarget, multi->sortedList-delta, multi->used, delta, insertpos, findCompare, TRUE);
436 if(insertpos > 0)
437 return( -1 );
438 insertpos = -insertpos - delta;
439
440 /* Check if the candidate is worse than the worst of the list */
441 if(((insertpos >= multi->size) && (multi->freeList[0] == 0)) ||
442 ((insertpos == multi->used) && (!allowSortedExpand ||
443 (multi->step_last >= multi->epszero))))
444 return( -1 );
445
446 #ifdef Paranoia
447 /* Do validation */
448 if((insertpos < 0) || (insertpos > multi->used))
449 return( -1 );
450 #endif
451
452 /* Define the target for storing the candidate;
453 Case 1: List is full and we must discard the previously worst candidate
454 Case 2: List is not full and we simply use the next free position */
455 if(multi->freeList[0] == 0)
456 targetrec = (pricerec *) multi->sortedList[multi->used-1].pvoidreal.ptr;
457 else {
458 delta = multi->freeList[0]--;
459 delta = multi->freeList[delta];
460 targetrec = &(multi->items[delta]);
461 }
462 }
463 else {
464 delta = multi->freeList[0]--;
465 delta = multi->freeList[delta];
466 targetrec = &(multi->items[delta]);
467 insertpos = multi->used;
468 }
469
470 /* Insert the new candidate record in the data store */
471 MEMCOPY(targetrec, candidate, 1);
472
473 /* Store the pointer data and handle tree cases:
474 Case 1: The list is unsorted and not full; simply append pointer to list,
475 Case 2: The list is sorted and full; insert the pointer by discarding previous last,
476 Case 3: The list is sorted and not full; shift the inferior items down, and increment count */
477 if((multi->used < multi->size) && (insertpos >= multi->used)) {
478 QS_append(multi->sortedList, insertpos, targetrec);
479 multi->used++;
480 }
481 else {
482 if(multi->used == multi->size)
483 QS_insert(multi->sortedList, insertpos, targetrec, multi->size-1); /* Discard previous last */
484 else {
485 QS_insert(multi->sortedList, insertpos, targetrec, multi->used); /* Keep previous last */
486 multi->used++;
487 }
488 }
489 multi->active = insertpos;
490
491 #ifdef Paranoia
492 if((insertpos >= multi->size) || (insertpos >= multi->used))
493 report(multi->lp, SEVERE, "addCandidateVar: Insertion point beyond limit!\n");
494 #endif
495
496 return( insertpos );
497 }
498
findImprovementVar(pricerec * current,pricerec * candidate,MYBOOL collectMP,int * candidatecount)499 STATIC MYBOOL findImprovementVar(pricerec *current, pricerec *candidate, MYBOOL collectMP, int *candidatecount)
500 /* PRIMAL: Find a variable to enter the basis
501 DUAL: Find a variable to leave the basis
502
503 Allowed variable set: Any pivot PRIMAL:larger or DUAL:smaller than threshold value of 0 */
504 {
505 MYBOOL Action = FALSE,
506 #ifdef ExtractedValidityTest
507 Accept = TRUE;
508 #else /* Check for validity and compare result with previous best */
509 Accept = validImprovementVar(candidate);
510 #endif
511 if(Accept) {
512 if(candidatecount != NULL)
513 (*candidatecount)++;
514 if(collectMP) {
515 if(addCandidateVar(candidate, current->lp->multivars, (findCompare_func *) compareImprovementQS, FALSE) < 0)
516 return(Action);
517 }
518 if(current->varno > 0)
519 Accept = (MYBOOL) (compareImprovementVar(current, candidate) > 0);
520 }
521
522 /* Apply candidate if accepted */
523 if(Accept) {
524 (*current) = *candidate;
525
526 /* Force immediate acceptance for Bland's rule using the primal simplex */
527 if(!candidate->isdual)
528 Action = (MYBOOL) (candidate->lp->_piv_rule_ == PRICER_FIRSTINDEX);
529 }
530 return(Action);
531 }
532
533 /* Bound flip variable accumulation routine */
collectMinorVar(pricerec * candidate,multirec * longsteps,MYBOOL isphase2,MYBOOL isbatch)534 STATIC MYBOOL collectMinorVar(pricerec *candidate, multirec *longsteps, MYBOOL isphase2, MYBOOL isbatch)
535 {
536 int inspos;
537
538 /* 1. Check for ratio and pivot validity (to have the extra flexibility that all
539 bound-flip candidates are also possible as basis-entering variables */
540 if(!validSubstitutionVar(candidate))
541 return( FALSE );
542
543 /* 2. If the free-list is empty we need to see if we have a better candidate,
544 and for this the candidate list has to be sorted by merit */
545 if(!isbatch &&
546 !longsteps->sorted && (longsteps->used > 1) &&
547 ((longsteps->freeList[0] == 0) ||
548 multi_truncatingvar(longsteps, candidate->varno) ||
549 (longsteps->step_last >= longsteps->epszero) )) {
550 longsteps->sorted = QS_execute(longsteps->sortedList, longsteps->used,
551 (findCompare_func *) compareSubstitutionQS, &inspos);
552 longsteps->dirty = (MYBOOL) (inspos > 0);
553 if(longsteps->dirty)
554 multi_recompute(longsteps, 0, isphase2, TRUE);
555 }
556
557 /* 3. Now handle three cases...
558 - Add to the list when the list is not full and there is opportunity for improvement,
559 - Check if we should replace an incumbent when the list is full,
560 - Check if we should replace an incumbent when the list is not full, there is no room
561 for improvement, but the current candidate is better than an incumbent. */
562 inspos = addCandidateVar(candidate, longsteps, (findCompare_func *) compareSubstitutionQS, TRUE);
563
564 /* 4. Recompute steps and objective, and (if relevant) determine if we
565 may be suboptimal in relation to an incumbent MILP solution. */
566 return( (MYBOOL) (inspos >= 0) &&
567 ((isbatch == TRUE) || multi_recompute(longsteps, inspos, isphase2, TRUE)) );
568 }
569
findSubstitutionVar(pricerec * current,pricerec * candidate,int * candidatecount)570 STATIC MYBOOL findSubstitutionVar(pricerec *current, pricerec *candidate, int *candidatecount)
571 /* PRIMAL: Find a variable to leave the basis
572 DUAL: Find a variable to enter the basis
573
574 Allowed variable set: "Equal-valued" smallest thetas! */
575 {
576 MYBOOL Action = FALSE,
577 #ifdef ExtractedValidityTest
578 Accept = TRUE;
579 #else /* Check for validity and comparison result with previous best */
580 Accept = validSubstitutionVar(candidate);
581 #endif
582 if(Accept) {
583 if(candidatecount != NULL)
584 (*candidatecount)++;
585 if(current->varno != 0)
586 Accept = (MYBOOL) (compareSubstitutionVar(current, candidate) > 0);
587 }
588
589 /* Apply candidate if accepted */
590 if(Accept) {
591 (*current) = *candidate;
592
593 /* Force immediate acceptance for Bland's rule using the dual simplex */
594 #ifdef ForceEarlyBlandRule
595 if(candidate->isdual)
596 Action = (MYBOOL) (candidate->lp->_piv_rule_ == PRICER_FIRSTINDEX);
597 #endif
598 }
599 return(Action);
600 }
601
602 /* Partial pricing management routines */
partial_createBlocks(lprec * lp,MYBOOL isrow)603 STATIC partialrec *partial_createBlocks(lprec *lp, MYBOOL isrow)
604 {
605 partialrec *blockdata;
606
607 blockdata = (partialrec *) calloc(1, sizeof(*blockdata));
608 blockdata->lp = lp;
609 blockdata->blockcount = 1;
610 blockdata->blocknow = 1;
611 blockdata->isrow = isrow;
612
613 return(blockdata);
614 }
partial_countBlocks(lprec * lp,MYBOOL isrow)615 STATIC int partial_countBlocks(lprec *lp, MYBOOL isrow)
616 {
617 partialrec *blockdata = IF(isrow, lp->rowblocks, lp->colblocks);
618
619 if(blockdata == NULL)
620 return( 1 );
621 else
622 return( blockdata->blockcount );
623 }
partial_activeBlocks(lprec * lp,MYBOOL isrow)624 STATIC int partial_activeBlocks(lprec *lp, MYBOOL isrow)
625 {
626 partialrec *blockdata = IF(isrow, lp->rowblocks, lp->colblocks);
627
628 if(blockdata == NULL)
629 return( 1 );
630 else
631 return( blockdata->blocknow );
632 }
partial_freeBlocks(partialrec ** blockdata)633 STATIC void partial_freeBlocks(partialrec **blockdata)
634 {
635 if((blockdata == NULL) || (*blockdata == NULL))
636 return;
637 FREE((*blockdata)->blockend);
638 FREE((*blockdata)->blockpos);
639 FREE(*blockdata);
640 }
641
642
643 /* Function to provide for left-right or right-left scanning of entering/leaving
644 variables; note that *end must have been initialized by the calling routine! */
makePriceLoop(lprec * lp,int * start,int * end,int * delta)645 STATIC void makePriceLoop(lprec *lp, int *start, int *end, int *delta)
646 {
647 int offset = is_piv_mode(lp, PRICE_LOOPLEFT);
648
649 if((offset) ||
650 (((lp->total_iter+offset) % 2 == 0) && is_piv_mode(lp, PRICE_LOOPALTERNATE))) {
651 *delta = -1; /* Step backwards - "left" */
652 swapINT(start, end);
653 lp->_piv_left_ = TRUE;
654 }
655 else {
656 *delta = 1; /* Step forwards - "right" */
657 lp->_piv_left_ = FALSE;
658 }
659 }
660
661 /* Routine to verify accuracy of the current basis factorization */
serious_facterror(lprec * lp,REAL * bvector,int maxcols,REAL tolerance)662 STATIC MYBOOL serious_facterror(lprec *lp, REAL *bvector, int maxcols, REAL tolerance)
663 {
664 int i, j, ib, ie, nz, nc;
665 REAL sum, tsum = 0, err = 0;
666 MATrec *mat = lp->matA;
667
668 if(bvector == 0)
669 bvector = lp->bsolveVal;
670 nc =0;
671 nz = 0;
672 for(i = 1; (i <= lp->rows) && (nc <= maxcols); i++) {
673
674 /* Do we have a non-slack variable? (we choose to skip slacks,
675 since they have "natural" good accuracy properties) */
676 j = lp->var_basic[i] - lp->rows;
677 if(j <= 0)
678 continue;
679 nc++;
680
681 /* Compute cross product for basic, non-slack column */
682 ib = mat->col_end[j-1];
683 ie = mat->col_end[j];
684 nz += ie - ib;
685 sum = get_OF_active(lp, j+lp->rows, bvector[0]);
686 for(; ib < ie; ib++)
687 sum += COL_MAT_VALUE(ib)*bvector[COL_MAT_ROWNR(ib)];
688
689 /* Catch high precision early, so we don't to uneccessary work */
690 tsum += sum;
691 SETMAX(err, fabs(sum));
692 if((tsum / nc > tolerance / 100) && (err < tolerance / 100))
693 break;
694 }
695 err /= mat->infnorm;
696 return( (MYBOOL) (err >= tolerance) );
697 }
698
699 /* Computation of reduced costs */
update_reducedcosts(lprec * lp,MYBOOL isdual,int leave_nr,int enter_nr,REAL * prow,REAL * drow)700 STATIC void update_reducedcosts(lprec *lp, MYBOOL isdual, int leave_nr, int enter_nr, REAL *prow, REAL *drow)
701 {
702 /* "Fast" update of the dual reduced cost vector; note that it must be called
703 after the pivot operation and only applies to a major "true" iteration */
704 int i;
705 REAL hold;
706
707 if(isdual) {
708 hold = -drow[enter_nr]/prow[enter_nr];
709 for(i=1; i <= lp->sum; i++)
710 if(!lp->is_basic[i]) {
711 if(i == leave_nr)
712 drow[i] = hold;
713 else {
714 drow[i] += hold*prow[i];
715 my_roundzero(drow[i], lp->epsmachine);
716 }
717 }
718 }
719 else
720 report(lp, SEVERE, "update_reducedcosts: Cannot update primal reduced costs!\n");
721 }
722
723
compute_reducedcosts(lprec * lp,MYBOOL isdual,int row_nr,int * coltarget,MYBOOL dosolve,REAL * prow,int * nzprow,REAL * drow,int * nzdrow,int roundmode)724 STATIC void compute_reducedcosts(lprec *lp, MYBOOL isdual, int row_nr, int *coltarget, MYBOOL dosolve,
725 REAL *prow, int *nzprow,
726 REAL *drow, int *nzdrow,
727 int roundmode)
728 {
729 REAL epsvalue = lp->epsvalue; /* Any larger value can produce a suboptimal result */
730 roundmode |= MAT_ROUNDRC;
731
732 if(isdual) {
733 bsolve_xA2(lp, coltarget,
734 row_nr, prow, epsvalue, nzprow, /* Calculate net sensitivity given a leaving variable */
735 0, drow, epsvalue, nzdrow, /* Calculate the net objective function values */
736 roundmode);
737 }
738 else {
739 REAL *bVector;
740
741 #if 1 /* Legacy mode, that is possibly a little faster */
742 if((lp->multivars == NULL) && (lp->P1extraDim == 0))
743 bVector = drow;
744 else
745 #endif
746 bVector = lp->bsolveVal;
747 if(dosolve) {
748 bsolve(lp, 0, bVector, lp->bsolveIdx, epsvalue*DOUBLEROUND, 1.0);
749 if(!isdual && (row_nr == 0) && (lp->improve & IMPROVE_SOLUTION) && !refactRecent(lp) &&
750 serious_facterror(lp, bVector, lp->rows, lp->epsvalue))
751 set_action(&lp->spx_action, ACTION_REINVERT);
752 }
753 prod_xA(lp, coltarget,
754 bVector, lp->bsolveIdx, epsvalue, 1.0,
755 drow, nzdrow, roundmode);
756 }
757 }
758
759
760 /* Primal: Prevent acceptance of an entering variable when the magnitude of
761 other candidates is also very small.
762 Dual: Prevent acceptance of a leaving variable when the magnitude of
763 other candidates is also very small.
764
765 Both of these cases are associated with numerical stalling, which we could
766 argue should be detected and handled by the stalling monitor routine. */
verify_stability(lprec * lp,MYBOOL isprimal,REAL xfeas,REAL sfeas,int nfeas)767 STATIC MYBOOL verify_stability(lprec *lp, MYBOOL isprimal, REAL xfeas, REAL sfeas, int nfeas)
768 {
769 MYBOOL testOK = TRUE;
770 return( testOK );
771
772 #if 1
773 /* Try to make dual feasibility as tight as possible */
774 if(!isprimal)
775 /* if(lp->P1extraVal == 0) */
776 {
777 xfeas /= (1+lp->rhsmax);
778 sfeas /= (1+lp->rhsmax);
779 }
780 #endif
781 xfeas = fabs(xfeas); /* Maximum (positive) infeasibility */
782 /* if(xfeas < lp->epspivot) { */
783 if(xfeas < lp->epssolution) {
784 REAL f;
785 sfeas = fabs(sfeas); /* Make sum of infeasibilities positive */
786 xfeas = (sfeas-xfeas)/nfeas; /* Average "residual" feasibility */
787 f = 1 + log10((REAL) nfeas); /* Some numerical complexity scalar */
788 /* Numerical errors can interact to cause non-convergence, and the
789 idea is to relax the tolerance to account for this and only
790 marginally weakening the (user-specified) tolerance. */
791 if((sfeas-xfeas) < f*lp->epsprimal)
792 testOK = FALSE;
793 }
794 return( testOK );
795 }
796
797
798 /* Find an entering column for the case that the specified basic variable
799 is fixed or zero - typically used for artificial variable elimination */
find_rowReplacement(lprec * lp,int rownr,REAL * prow,int * nzprow)800 STATIC int find_rowReplacement(lprec *lp, int rownr, REAL *prow, int *nzprow)
801 /* The logic in this section generally follows Chvatal: Linear Programming, p. 130
802 Basically, the function is a specialized coldual(). */
803 {
804 int i, bestindex;
805 REAL bestvalue;
806
807 /* Solve for "local reduced cost" */
808 set_action(&lp->piv_strategy, PRICE_FORCEFULL);
809 compute_reducedcosts(lp, TRUE, rownr, NULL, TRUE,
810 prow, nzprow, NULL, NULL, MAT_ROUNDDEFAULT);
811 clear_action(&lp->piv_strategy, PRICE_FORCEFULL);
812
813 /* Find a suitably non-singular variable to enter ("most orthogonal") */
814 bestindex = 0;
815 bestvalue = 0;
816 for(i = 1; i <= lp->sum-abs(lp->P1extraDim); i++) {
817 if(!lp->is_basic[i] && !is_fixedvar(lp, i) &&
818 (fabs(prow[i]) > bestvalue)) {
819 bestindex = i;
820 bestvalue = fabs(prow[i]);
821 }
822 }
823
824 /* Prepare to update inverse and pivot/iterate (compute Bw=a) */
825 if(i > lp->sum-abs(lp->P1extraDim))
826 bestindex = 0;
827 else
828 fsolve(lp, bestindex, prow, nzprow, lp->epsmachine, 1.0, TRUE);
829
830 return( bestindex );
831 }
832
833 /* Find the primal simplex entering non-basic column variable */
colprim(lprec * lp,REAL * drow,int * nzdrow,MYBOOL skipupdate,int partialloop,int * candidatecount,MYBOOL updateinfeas,REAL * xviol)834 STATIC int colprim(lprec *lp, REAL *drow, int *nzdrow, MYBOOL skipupdate, int partialloop, int *candidatecount, MYBOOL updateinfeas, REAL *xviol)
835 {
836 int i, ix, iy, iz, ninfeas, nloop = 0;
837 REAL f, sinfeas, xinfeas, epsvalue = lp->epsdual;
838 pricerec current, candidate;
839 MYBOOL collectMP = FALSE;
840 int *coltarget = NULL;
841
842 /* Identify pivot column according to pricing strategy; set
843 entering variable initial threshold reduced cost value to "0" */
844 current.pivot = lp->epsprimal; /* Minimum acceptable improvement */
845 current.varno = 0;
846 current.lp = lp;
847 current.isdual = FALSE;
848 candidate.lp = lp;
849 candidate.isdual = FALSE;
850 *candidatecount = 0;
851
852 /* Update local value of pivot setting and determine active multiple pricing set */
853 lp->_piv_rule_ = get_piv_rule(lp);
854 doLoop:
855 nloop++;
856 if((lp->multivars != NULL) && ((lp->simplex_mode & SIMPLEX_PRIMAL_PRIMAL) != 0)) {
857 collectMP = multi_mustupdate(lp->multivars);
858 if(collectMP) {
859 multi_restart(lp->multivars);
860 coltarget = NULL;
861 }
862 else
863 coltarget = multi_indexSet(lp->multivars, FALSE);
864 }
865
866 /* Compute reduced costs c - c*Inv(B), if necessary
867 (i.e. the previous iteration was not a "minor" iteration/bound flip) */
868 if(!skipupdate) {
869 #ifdef UsePrimalReducedCostUpdate
870 /* Recompute from scratch only at the beginning, otherwise update */
871 if((lp->current_iter > 0) && (refactRecent(lp) == AUTOMATIC))
872 #endif
873 compute_reducedcosts(lp, FALSE, 0, coltarget, (MYBOOL) ((nloop <= 1) || (partialloop > 1)),
874 NULL, NULL,
875 drow, nzdrow,
876 MAT_ROUNDDEFAULT);
877 }
878
879 /* Loop over active partial column set; we presume that reduced costs
880 have only been updated for columns in the active partial range. */
881 ix = 1;
882 iy = nzdrow[0];
883 ninfeas = 0;
884 xinfeas = 0;
885 sinfeas = 0;
886 makePriceLoop(lp, &ix, &iy, &iz);
887 iy *= iz;
888 for(; ix*iz <= iy; ix += iz) {
889 i = nzdrow[ix];
890 #if 0 /* Not necessary since we masked them out in compute_reducedcosts() */
891 if(i > lp->sum-abs(lp->P1extraDim))
892 continue;
893 #endif
894
895 /* Check if the pivot candidate is on the block-list */
896 if(lp->rejectpivot[0] > 0) {
897 int kk;
898 for(kk = 1; (kk <= lp->rejectpivot[0]) && (i != lp->rejectpivot[kk]); kk++);
899 if(kk <= lp->rejectpivot[0])
900 continue;
901 }
902
903 /* Retrieve the applicable reduced cost - threshold should not be smaller than 0 */
904 f = my_chsign(lp->is_lower[i], drow[i]);
905 if(f <= epsvalue)
906 continue;
907
908 /* Find entering variable according to strategy (largest positive f) */
909 ninfeas++;
910 SETMAX(xinfeas, f);
911 sinfeas += f;
912 candidate.pivot = normalizeEdge(lp, i, f, FALSE);
913 candidate.varno = i;
914 if(findImprovementVar(¤t, &candidate, collectMP, candidatecount))
915 break;
916 }
917
918 /* Check if we should loop again after a multiple pricing update */
919 if(lp->multivars != NULL) {
920 if(collectMP) {
921 if(!lp->multivars->sorted)
922 lp->multivars->sorted = QS_execute(lp->multivars->sortedList, lp->multivars->used,
923 (findCompare_func *) compareImprovementQS, NULL);
924 coltarget = multi_indexSet(lp->multivars, TRUE);
925 }
926 else if((current.varno == 0) && (lp->multivars->retries == 0)) {
927 ix = partial_blockStart(lp, FALSE);
928 iy = partial_blockEnd(lp, FALSE);
929 lp->multivars->used = 0;
930 lp->multivars->retries++;
931 goto doLoop;
932 }
933 /* Shrink the candidate list */
934 lp->multivars->retries = 0;
935 if(current.varno != 0)
936 multi_removevar(lp->multivars, current.varno);
937 }
938
939 /* Check for optimality */
940 if(xviol != NULL)
941 *xviol = xinfeas;
942 if(updateinfeas)
943 lp->suminfeas = fabs(sinfeas);
944 if((lp->multivars == NULL) && (current.varno > 0) &&
945 !verify_stability(lp, TRUE, xinfeas, sinfeas, ninfeas))
946 current.varno = 0;
947
948 /* Produce statistics */
949 if(lp->spx_trace) {
950 if(current.varno > 0)
951 report(lp, DETAILED, "colprim: Column %d reduced cost = " RESULTVALUEMASK "\n",
952 current.varno, current.pivot);
953 else
954 report(lp, DETAILED, "colprim: No positive reduced costs found, optimality!\n");
955 }
956
957 return( current.varno );
958 } /* colprim */
959
960 /* Find the primal simplex leaving basic column variable */
rowprim(lprec * lp,int colnr,LREAL * theta,REAL * pcol,int * nzpcol,MYBOOL forceoutEQ,REAL * xviol)961 STATIC int rowprim(lprec *lp, int colnr, LREAL *theta, REAL *pcol, int *nzpcol, MYBOOL forceoutEQ, REAL *xviol)
962 {
963 int i, ii, iy, iz, Hpass, k, *nzlist;
964 LREAL f, savef;
965 REAL Heps, Htheta, Hlimit, epsvalue, epspivot, p;
966 pricerec current, candidate;
967 MYBOOL isupper = !lp->is_lower[colnr], HarrisTwoPass = FALSE;
968
969 p = 0.0;
970
971 /* Update local value of pivot setting */
972 lp->_piv_rule_ = get_piv_rule(lp);
973 if(nzpcol == NULL)
974 nzlist = (int *) mempool_obtainVector(lp->workarrays, lp->rows+1, sizeof(*nzlist));
975 else
976 nzlist = nzpcol;
977
978 /* Find unconditional non-zeros and optionally compute relative size of epspivot */
979 epspivot = lp->epspivot;
980 epsvalue = lp->epsvalue;
981 Hlimit = 0;
982 Htheta = 0;
983 k = 0;
984 for(i = 1; i <= lp->rows; i++) {
985 p = fabs(pcol[i]);
986 if(p > Hlimit)
987 Hlimit = p;
988 if(p > epsvalue) {
989 k++;
990 nzlist[k] = i;
991 SETMAX(Htheta, p);
992 }
993 #ifdef Paranoia
994 else {
995 if(lp->spx_trace)
996 report(lp, FULL, "rowprim: Row %d with pivot " RESULTVALUEMASK " rejected as too small\n",
997 i, p);
998 }
999 #endif
1000 }
1001 if(xviol != NULL)
1002 *xviol = Htheta;
1003 Htheta = 0;
1004
1005 /* Update non-zero list based on the new pivot threshold */
1006 #ifdef UseRelativePivot_Primal
1007 /* epspivot *= sqrt(lp->matA->dynrange) / lp->matA->infnorm; */
1008 epspivot /= MAX(1, sqrt(lp->matA->colmax[colnr]));
1009 iy = k;
1010 k = 0;
1011 p = 0;
1012 for(ii = 1; ii <= iy; ii++) {
1013 i = nzlist[ii];
1014 p = fabs(pcol[i]);
1015
1016 /* Compress the list of valid alternatives, if appropriate */
1017 if(p > epspivot) {
1018 k++;
1019 nzlist[k] = i;
1020 }
1021 #ifdef Paranoia
1022 else {
1023 if(lp->spx_trace)
1024 report(lp, FULL, "rowprim: Row %d with pivot " RESULTVALUEMASK " rejected as too small\n",
1025 i, p);
1026 }
1027 #endif
1028 }
1029 #endif
1030
1031 /* Initialize counters */
1032 nzlist[0] = k;
1033 k = 0;
1034
1035 Retry:
1036 k++;
1037 HarrisTwoPass = is_piv_mode(lp, PRICE_HARRISTWOPASS);
1038 if(HarrisTwoPass)
1039 Hpass = 1;
1040 else
1041 Hpass = 2;
1042 current.theta = lp->infinite;
1043 current.pivot = 0;
1044 current.varno = 0;
1045 current.isdual = FALSE;
1046 current.epspivot = epspivot;
1047 current.lp = lp;
1048 candidate.epspivot = epspivot;
1049 candidate.isdual = FALSE;
1050 candidate.lp = lp;
1051 savef = 0;
1052 for(; Hpass <= 2; Hpass++) {
1053 Htheta = lp->infinite;
1054 if(Hpass == 1) {
1055 Hlimit = lp->infinite; /* Don't apply any limit in the first pass */
1056 Heps = epspivot/lp->epsprimal; /* Scaled to lp->epsprimal used in compute_theta() */
1057 }
1058 else {
1059 Hlimit = current.theta; /* This is the smallest Theta of the first pass */
1060 Heps = 0.0;
1061 }
1062 current.theta = lp->infinite;
1063 current.pivot = 0;
1064 current.varno = 0;
1065 savef = 0;
1066
1067 ii = 1;
1068 iy = nzlist[0];
1069 makePriceLoop(lp, &ii, &iy, &iz);
1070 iy *= iz;
1071 for(; ii*iz <= iy; ii += iz) {
1072 i = nzlist[ii];
1073 f = pcol[i];
1074 candidate.theta = f;
1075 candidate.pivot = f;
1076 candidate.varno = i;
1077
1078 /*i =*/ compute_theta(lp, i, &candidate.theta, isupper,
1079 my_if(lp->upbo[lp->var_basic[i]] < lp->epsprimal, Heps/10, Heps), TRUE);
1080
1081 if(fabs(candidate.theta) >= lp->infinite) {
1082 savef = f;
1083 candidate.theta = 2*lp->infinite;
1084 continue;
1085 }
1086
1087 /* Find the candidate leaving variable according to strategy (smallest theta) */
1088 if((Hpass == 2) && (candidate.theta > Hlimit))
1089 continue;
1090
1091 /* Give a slight preference to fixed variables (mainly equality slacks) */
1092 if(forceoutEQ) {
1093 p = candidate.pivot;
1094 if(lp->upbo[lp->var_basic[i]] < lp->epsprimal) {
1095 /* Give an extra early boost to equality slack elimination, if specified */
1096 if(forceoutEQ == AUTOMATIC)
1097 candidate.pivot *= 1.0+lp->epspivot;
1098 else
1099 candidate.pivot *= 10.0;
1100
1101 }
1102 }
1103 if(HarrisTwoPass) {
1104 f = candidate.theta;
1105 if(Hpass == 2)
1106 candidate.theta = 1;
1107 if(findSubstitutionVar(¤t, &candidate, NULL))
1108 break;
1109 if((Hpass == 2) && (current.varno == candidate.varno))
1110 Htheta = f;
1111 }
1112 else
1113 if(findSubstitutionVar(¤t, &candidate, NULL))
1114 break;
1115 /* Restore temporarily modified pivot */
1116 if(forceoutEQ && (current.varno == candidate.varno))
1117 current.pivot = p;
1118 }
1119 }
1120 if(HarrisTwoPass)
1121 current.theta = Htheta;
1122
1123 /* Handle case of no available leaving variable */
1124 if(current.varno == 0) {
1125 if(lp->upbo[colnr] >= lp->infinite) {
1126 /* Optionally try again with reduced pivot threshold level */
1127 if(k < 2) {
1128 epspivot = epspivot / 10;
1129 goto Retry;
1130 }
1131 }
1132 else {
1133 #if 1
1134 i = 1;
1135 while((pcol[i] >= 0) && (i <= lp->rows))
1136 i++;
1137 if(i > lp->rows) { /* Empty column with upper bound! */
1138 lp->is_lower[colnr] = !lp->is_lower[colnr];
1139 /* lp->is_lower[colnr] = FALSE; */
1140 lp->rhs[0] += lp->upbo[colnr]*pcol[0];
1141 }
1142 else /* if(pcol[i]<0) */
1143 {
1144 current.varno = i;
1145 }
1146 #endif
1147 }
1148 }
1149 else if(current.theta >= lp->infinite) {
1150 report(lp, IMPORTANT, "rowprim: Numeric instability pcol[%d] = %g, rhs[%d] = %g, upbo = %g\n",
1151 current.varno, savef, current.varno, lp->rhs[current.varno],
1152 lp->upbo[lp->var_basic[current.varno]]);
1153 }
1154
1155 /* Return working array to pool */
1156 if(nzpcol == NULL)
1157 mempool_releaseVector(lp->workarrays, (char *) nzlist, FALSE);
1158
1159 if(lp->spx_trace)
1160 report(lp, DETAILED, "row_prim: %d, pivot size = " RESULTVALUEMASK "\n",
1161 current.varno, current.pivot);
1162
1163 /* *theta = current.theta; */
1164 *theta = fabs(current.theta);
1165
1166 return(current.varno);
1167 } /* rowprim */
1168
1169
1170 /* Find the dual simplex leaving basic variable */
rowdual(lprec * lp,REAL * rhvec,MYBOOL forceoutEQ,MYBOOL updateinfeas,REAL * xviol)1171 STATIC int rowdual(lprec *lp, REAL *rhvec, MYBOOL forceoutEQ, MYBOOL updateinfeas, REAL *xviol)
1172 {
1173 int k, i, iy, iz, ii, ninfeas;
1174 register REAL rh;
1175 REAL up, lo = 0,
1176 epsvalue, sinfeas, xinfeas;
1177 pricerec current, candidate;
1178 MYBOOL collectMP = FALSE;
1179
1180 /* Initialize */
1181 if(rhvec == NULL)
1182 rhvec = lp->rhs;
1183 epsvalue = lp->epsdual;
1184 current.pivot = -epsvalue; /* Initialize leaving variable threshold; "less than 0" */
1185 current.theta = 0;
1186 current.varno = 0;
1187 current.isdual = TRUE;
1188 current.lp = lp;
1189 candidate.isdual = TRUE;
1190 candidate.lp = lp;
1191
1192 /* Loop over active partial row set */
1193 if(is_action(lp->piv_strategy, PRICE_FORCEFULL)) {
1194 k = 1;
1195 iy = lp->rows;
1196 }
1197 else {
1198 k = partial_blockStart(lp, TRUE);
1199 iy = partial_blockEnd(lp, TRUE);
1200 }
1201 ninfeas = 0;
1202 xinfeas = 0;
1203 sinfeas = 0;
1204 makePriceLoop(lp, &k, &iy, &iz);
1205 iy *= iz;
1206 for(; k*iz <= iy; k += iz) {
1207
1208 /* Map loop variable to target */
1209 i = k;
1210
1211 /* Check if the pivot candidate is on the block-list */
1212 if(lp->rejectpivot[0] > 0) {
1213 int kk;
1214 for(kk = 1; (kk <= lp->rejectpivot[0]) && (i != lp->rejectpivot[kk]); kk++);
1215 if(kk <= lp->rejectpivot[0])
1216 continue;
1217 }
1218
1219 /* Set local variables - express violation as a negative number */
1220 ii = lp->var_basic[i];
1221 up = lp->upbo[ii];
1222 lo = 0;
1223 rh = rhvec[i];
1224 if(rh > up)
1225 rh = up - rh;
1226 else
1227 rh -= lo;
1228 up -= lo;
1229
1230 /* Analyze relevant constraints ...
1231 KE version skips uninteresting alternatives and gives a noticeable speedup */
1232 /* if((rh < -epsvalue*sqrt(lp->matA->rowmax[i])) || */
1233 if((rh < -epsvalue) ||
1234 ((forceoutEQ == TRUE) && (up < epsvalue))) { /* It causes instability to remove the "TRUE" test */
1235
1236 /* Accumulate stats */
1237 ninfeas++;
1238 SETMIN(xinfeas, rh);
1239 sinfeas += rh;
1240
1241 /* Give a slight preference to fixed variables (mainly equality slacks) */
1242 if(up < epsvalue) {
1243 /* Break out immediately if we are directed to force slacks out of the basis */
1244 if(forceoutEQ == TRUE) {
1245 current.varno = i;
1246 current.pivot = -1;
1247 break;
1248 }
1249 /* Give an extra early boost to equality slack elimination, if specified */
1250 if(forceoutEQ == AUTOMATIC)
1251 rh *= 10.0;
1252 else /* .. or just the normal. marginal boost */
1253 rh *= 1.0+lp->epspivot;
1254 }
1255
1256 /* Select leaving variable according to strategy (the most negative/largest violation) */
1257 candidate.pivot = normalizeEdge(lp, i, rh, TRUE);
1258 candidate.varno = i;
1259 if(findImprovementVar(¤t, &candidate, collectMP, NULL))
1260 break;
1261 }
1262 }
1263
1264 /* Verify infeasibility */
1265 if(updateinfeas)
1266 lp->suminfeas = fabs(sinfeas);
1267 if((ninfeas > 1) &&
1268 !verify_stability(lp, FALSE, xinfeas, sinfeas, ninfeas)) {
1269 report(lp, IMPORTANT, "rowdual: Check for reduced accuracy and tolerance settings.\n");
1270 current.varno = 0;
1271 }
1272
1273 /* Produce statistics */
1274 if(lp->spx_trace) {
1275 report(lp, NORMAL, "rowdual: Infeasibility sum " RESULTVALUEMASK " in %7d constraints.\n",
1276 sinfeas, ninfeas);
1277 if(current.varno > 0) {
1278 report(lp, DETAILED, "rowdual: rhs[%d] = " RESULTVALUEMASK "\n",
1279 current.varno, lp->rhs[current.varno]);
1280 }
1281 else
1282 report(lp, FULL, "rowdual: Optimality - No primal infeasibilities found\n");
1283 }
1284 if(xviol != NULL)
1285 *xviol = fabs(xinfeas);
1286
1287 return(current.varno);
1288 } /* rowdual */
1289
1290
longdual_testset(lprec * lp,int which,int rownr,REAL * prow,int * nzprow,REAL * drow,int * nzdrow)1291 STATIC void longdual_testset(lprec *lp, int which, int rownr, REAL *prow, int *nzprow,
1292 REAL *drow, int *nzdrow)
1293 {
1294 int i,j;
1295 REAL F = lp->infinite;
1296 if(which == 0) { /* Maros Example-1 - raw data */
1297 j = 1; i = lp->rows+j; lp->upbo[i] = 0; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = 2; drow[i] = -1;
1298 j = 2; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = -2; drow[i] = 2;
1299 j = 3; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = 1; drow[i] = 5;
1300 j = 4; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = 3; drow[i] = -6;
1301 j = 5; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = -4; drow[i] = -2;
1302 j = 6; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = -1; drow[i] = 0;
1303 j = 7; i = lp->rows+j; lp->upbo[i] = 2; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = 1; drow[i] = 0;
1304 j = 8; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = -2; drow[i] = 0;
1305 j = 9; i = lp->rows+j; lp->upbo[i] = 5; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = -1; drow[i] = 4;
1306 j = 10; i = lp->rows+j; lp->upbo[i] = F; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = -2; drow[i] = 10;
1307 nzprow[0] = i-lp->rows;
1308 lp->rhs[rownr] = -11;
1309 lp->upbo[lp->var_basic[rownr]] = F;
1310 lp->rhs[0] = 1;
1311 }
1312 else if(which == 1) { /* Maros Example-1 - presorted in correct order */
1313 j = 1; i = lp->rows+j; lp->upbo[i] = 0; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = 2; drow[i] = -1;
1314 j = 2; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = 1; drow[i] = 5;
1315 j = 3; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = -4; drow[i] = -2;
1316 j = 4; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = -2; drow[i] = 0;
1317
1318 j = 5; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = -1; drow[i] = 0;
1319 j = 6; i = lp->rows+j; lp->upbo[i] = 2; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = 1; drow[i] = 0;
1320 j = 7; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = -2; drow[i] = 2;
1321 j = 8; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = 3; drow[i] = -6;
1322 j = 9; i = lp->rows+j; lp->upbo[i] = 5; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = -1; drow[i] = 4;
1323 j = 10; i = lp->rows+j; lp->upbo[i] = F; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = -2; drow[i] = 10;
1324 nzprow[0] = i-lp->rows;
1325 lp->rhs[rownr] = -11;
1326 lp->upbo[lp->var_basic[rownr]] = F;
1327 lp->rhs[0] = 1;
1328 }
1329
1330 else if(which == 10) { /* Maros Example-2 - raw data */
1331 j = 1; i = lp->rows+j; lp->upbo[i] = 5; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = -2; drow[i] = 2;
1332 j = 2; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = 3; drow[i] = 3;
1333 j = 3; i = lp->rows+j; lp->upbo[i] = 1; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = -2; drow[i] = 0;
1334 j = 4; i = lp->rows+j; lp->upbo[i] = 2; lp->is_lower[i] = FALSE; nzprow[j] = i; prow[i] = -1; drow[i] = -2;
1335 j = 5; i = lp->rows+j; lp->upbo[i] = 2; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = 1; drow[i] = 0;
1336 j = 6; i = lp->rows+j; lp->upbo[i] = F; lp->is_lower[i] = TRUE; nzprow[j] = i; prow[i] = 3; drow[i] = 9;
1337 nzprow[0] = i-lp->rows;
1338 lp->rhs[rownr] = 14;
1339 lp->upbo[lp->var_basic[rownr]] = 2;
1340 lp->rhs[0] = 6;
1341 }
1342 }
1343
1344
1345 /* Find the dual simplex entering non-basic variable */
coldual(lprec * lp,int row_nr,REAL * prow,int * nzprow,REAL * drow,int * nzdrow,MYBOOL dualphase1,MYBOOL skipupdate,int * candidatecount,REAL * xviol)1346 STATIC int coldual(lprec *lp, int row_nr, REAL *prow, int *nzprow,
1347 REAL *drow, int *nzdrow,
1348 MYBOOL dualphase1, MYBOOL skipupdate,
1349 int *candidatecount, REAL *xviol)
1350 {
1351 int i, iy, iz, ix, k, nbound;
1352 LREAL w, g, quot;
1353 REAL viol, p, epspivot = lp->epspivot;
1354 #ifdef MachinePrecRoundRHS
1355 REAL epsvalue = lp->epsmachine;
1356 #else
1357 REAL epsvalue = lp->epsvalue;
1358 #endif
1359 pricerec current, candidate;
1360 MYBOOL isbatch = FALSE, /* Requires that lp->longsteps->size > lp->sum */
1361 dolongsteps = (MYBOOL) (lp->longsteps != NULL);
1362
1363 /* Initialize */
1364 if(dolongsteps && !dualphase1)
1365 dolongsteps = AUTOMATIC; /* Sets Phase1 = TRUE, Phase2 = AUTOMATIC */
1366 current.theta = lp->infinite;
1367 current.pivot = 0;
1368 current.varno = 0;
1369 current.epspivot = epspivot;
1370 current.isdual = TRUE;
1371 current.lp = lp;
1372 candidate.epspivot = epspivot;
1373 candidate.isdual = TRUE;
1374 candidate.lp = lp;
1375 *candidatecount = 0;
1376
1377 /* Compute reduced costs */
1378 if(!skipupdate) {
1379 #ifdef UseDualReducedCostUpdate
1380 /* Recompute from scratch only at the beginning, otherwise update */
1381 if((lp->current_iter > 0) && (refactRecent(lp) < AUTOMATIC))
1382 compute_reducedcosts(lp, TRUE, row_nr, NULL, TRUE,
1383 prow, nzprow,
1384 NULL, NULL,
1385 MAT_ROUNDDEFAULT);
1386 else
1387 #endif
1388 compute_reducedcosts(lp, TRUE, row_nr, NULL, TRUE,
1389 prow, nzprow,
1390 drow, nzdrow,
1391 MAT_ROUNDDEFAULT);
1392 }
1393
1394 #if 0
1395 /* Override all above to do in-line testing with fixed test set */
1396 if(lp->rows > 1 && lp->columns > 10)
1397 longdual_testset(lp, 10, row_nr, prow, nzprow, drow, nzdrow);
1398 #endif
1399
1400 /* Compute the current violation of the bounds of the outgoing variable,
1401 negative for violation of lower bound, positive for upper bound violation.
1402 (Basic variables are always lower-bounded, by lp_solve convention) */
1403 g = 1;
1404 viol = lp->rhs[row_nr];
1405 if(viol > 0) { /* Check if the leaving variable is >= its upper bound */
1406 p = lp->upbo[lp->var_basic[row_nr]];
1407 if(p < lp->infinite) {
1408 viol -= p;
1409 my_roundzero(viol, epsvalue);
1410 if(viol > 0)
1411 g = -1;
1412 }
1413 /* Do validation of numerics */
1414 if(g == 1) {
1415 if(viol >= lp->infinite) {
1416 report(lp, IMPORTANT, "coldual: Large basic solution value %g at iter %.0f indicates numerical instability\n",
1417 lp->rhs[row_nr], (double) get_total_iter(lp));
1418 lp->spx_status = NUMFAILURE;
1419 return( 0 );
1420
1421 }
1422 if(skipupdate)
1423 report(lp, DETAILED, "coldual: Inaccurate bound-flip accuracy at iter %.0f\n",
1424 (double) get_total_iter(lp));
1425 else
1426 report(lp, SEVERE, "coldual: Leaving variable %d does not violate bounds at iter %.0f\n",
1427 row_nr, (double) get_total_iter(lp));
1428 return( -1 );
1429 }
1430 }
1431
1432 /* Update local value of pivot setting */
1433 lp->_piv_rule_ = get_piv_rule(lp);
1434
1435 /* Condense list of relevant targets */
1436 p = 0;
1437 k = 0;
1438 nbound = 0;
1439 ix = 1;
1440 iy = nzprow[0];
1441 for(ix = 1; ix <= iy; ix++) {
1442 i = nzprow[ix];
1443 w = prow[i] * g; /* Change sign if upper bound of the leaving variable is violated */
1444 w *= 2*lp->is_lower[i] - 1; /* Change sign if the non-basic variable is currently upper-bounded */
1445
1446 /* Check if the candidate is worth using for anything */
1447 if(w < -epsvalue) {
1448 /* Tally bounded variables */
1449 if(lp->upbo[i] < lp->infinite)
1450 nbound++;
1451
1452 /* Update the nz-index */
1453 k++;
1454 nzprow[k] = i;
1455 SETMAX(p, -w);
1456 }
1457 #ifdef Paranoia
1458 else {
1459 if(lp->spx_trace) {
1460 report(lp, FULL, "coldual: Candidate variable prow[%d] rejected with %g too small\n",
1461 i, w);
1462 }
1463 }
1464 #endif
1465
1466 }
1467 nzprow[0] = k;
1468 if(xviol != NULL)
1469 *xviol = p;
1470
1471 #ifdef UseRelativePivot_Dual
1472 /* epspivot *= sqrt(lp->matA->dynrange) / lp->matA->infnorm; */
1473 epspivot /= MAX(1, sqrt(lp->matA->rowmax[row_nr]));
1474 #endif
1475 current.epspivot = epspivot;
1476 candidate.epspivot = epspivot;
1477
1478 /* Initialize the long-step structures if indicated */
1479 if(dolongsteps) {
1480 if((nzprow[0] <= 1) || (nbound == 0)) { /* Don't bother */
1481 dolongsteps = FALSE;
1482 lp->longsteps->indexSet[0] = 0;
1483 }
1484 else {
1485 multi_restart(lp->longsteps);
1486 multi_valueInit(lp->longsteps, g*viol, lp->rhs[0]);
1487 }
1488 }
1489
1490 /* Loop over all entering column candidates */
1491 ix = 1;
1492 iy = nzprow[0];
1493 makePriceLoop(lp, &ix, &iy, &iz);
1494 iy *= iz;
1495 for(; ix*iz <= iy; ix += iz) {
1496 i = nzprow[ix];
1497
1498 /* Compute the dual ratio (prow = w and drow = cbar in Chvatal's "nomenclatura") */
1499 w = prow[i] * g; /* Change sign if upper bound of the leaving variable is violated */
1500 quot = -drow[i] / w; /* Remember this sign-reversal in multi_recompute! */
1501
1502 /* Apply the selected pivot strategy (smallest theta) */
1503 candidate.theta = quot; /* Note that abs() is applied in findSubstitutionVar */
1504 candidate.pivot = w;
1505 candidate.varno = i;
1506
1507 /* Collect candidates for minor iterations/bound flips */
1508 if(dolongsteps) {
1509 if(isbatch && (ix == iy))
1510 isbatch = AUTOMATIC;
1511 if(collectMinorVar(&candidate, lp->longsteps, (MYBOOL) (dolongsteps == AUTOMATIC), isbatch) &&
1512 lp->spx_trace)
1513 report(lp, DETAILED, "coldual: Long-dual break point with %d bound-flip variables\n",
1514 lp->longsteps->used);
1515 if(lp->spx_status == FATHOMED)
1516 return( 0 );
1517 }
1518
1519 /* We have a candidate for entering the basis; check if it is better than the incumbent */
1520 else if(findSubstitutionVar(¤t, &candidate, candidatecount))
1521 break;
1522 }
1523
1524 /* Set entering variable and long-step bound swap variables */
1525 if(dolongsteps) {
1526 *candidatecount = lp->longsteps->used;
1527 i = multi_enteringvar(lp->longsteps, NULL, 3);
1528 }
1529 else
1530 i = current.varno;
1531
1532 if(lp->spx_trace)
1533 report(lp, NORMAL, "coldual: Entering column %d, reduced cost %g, pivot value %g, bound swaps %d\n",
1534 i, drow[i], prow[i], multi_used(lp->longsteps));
1535
1536 return( i );
1537 } /* coldual */
1538
1539
normalizeEdge(lprec * lp,int item,REAL edge,MYBOOL isdual)1540 INLINE REAL normalizeEdge(lprec *lp, int item, REAL edge, MYBOOL isdual)
1541 {
1542 #if 1
1543 /* Don't use the pricer "close to home", since this can possibly
1544 worsen the final feasibility picture (mainly a Devex issue?) */
1545 if(fabs(edge) > lp->epssolution)
1546 #endif
1547 edge /= getPricer(lp, item, isdual);
1548 if((lp->piv_strategy & PRICE_RANDOMIZE) != 0)
1549 edge *= (1.0-PRICER_RANDFACT) + PRICER_RANDFACT*rand_uniform(lp, 1.0);
1550 return( edge );
1551
1552 }
1553
1554 /* Support routines for block detection and partial pricing */
partial_findBlocks(lprec * lp,MYBOOL autodefine,MYBOOL isrow)1555 STATIC int partial_findBlocks(lprec *lp, MYBOOL autodefine, MYBOOL isrow)
1556 {
1557 int i, jj, n, nb, ne, items;
1558 REAL hold, biggest, *sum = NULL;
1559 MATrec *mat = lp->matA;
1560 partialrec *blockdata;
1561
1562 if(!mat_validate(mat))
1563 return( 1 );
1564
1565 blockdata = IF(isrow, lp->rowblocks, lp->colblocks);
1566 items = IF(isrow, lp->rows, lp->columns);
1567 allocREAL(lp, &sum, items+1, FALSE);
1568
1569 /* Loop over items and compute the average column index for each */
1570 sum[0] = 0;
1571 for(i = 1; i <= items; i++) {
1572 n = 0;
1573 if(isrow) {
1574 nb = mat->row_end[i-1];
1575 ne = mat->row_end[i];
1576 }
1577 else {
1578 nb = mat->col_end[i-1];
1579 ne = mat->col_end[i];
1580 }
1581 n = ne-nb;
1582 sum[i] = 0;
1583 if(n > 0) {
1584 if(isrow)
1585 for(jj = nb; jj < ne; jj++)
1586 sum[i] += ROW_MAT_COLNR(jj);
1587 else
1588 for(jj = nb; jj < ne; jj++)
1589 sum[i] += COL_MAT_ROWNR(jj);
1590 sum[i] /= n;
1591 }
1592 else
1593 sum[i] = sum[i-1];
1594 }
1595
1596 /* Loop over items again, find largest difference and make monotone */
1597 hold = 0;
1598 biggest = 0;
1599 for(i = 2; i <= items; i++) {
1600 hold = sum[i] - sum[i-1];
1601 if(hold > 0) {
1602 if(hold > biggest)
1603 biggest = hold;
1604 }
1605 else
1606 hold = 0;
1607 sum[i-1] = hold;
1608 }
1609
1610 /* Loop over items again and find differences exceeding threshold;
1611 the discriminatory power of this routine depends strongly on the
1612 magnitude of the scaling factor - from empirical evidence > 0.9 */
1613 biggest = MAX(1, 0.9*biggest);
1614 n = 0;
1615 nb = 0;
1616 ne = 0;
1617 for(i = 1; i < items; i++)
1618 if(sum[i] > biggest) {
1619 ne += i-nb; /* Compute sum of index gaps between maxima */
1620 nb = i;
1621 n++; /* Increment count */
1622 }
1623
1624 /* Clean up */
1625 FREE(sum);
1626
1627 /* Require that the maxima are spread "nicely" across the columns,
1628 otherwise return that there is only one monolithic block.
1629 (This is probably an area for improvement in the logic!) */
1630 if(n > 0) {
1631 ne /= n; /* Average index gap between maxima */
1632 i = IF(isrow, lp->columns, lp->rows);
1633 nb = i / ne; /* Another estimated block count */
1634 if(abs(nb - n) > 2) /* Probably Ok to require equality (nb==n)*/
1635 n = 1;
1636 else if(autodefine) /* Generate row/column break-indeces for partial pricing */
1637 set_partialprice(lp, nb, NULL, isrow);
1638 }
1639 else
1640 n = 1;
1641
1642 return( n );
1643 }
partial_blockStart(lprec * lp,MYBOOL isrow)1644 STATIC int partial_blockStart(lprec *lp, MYBOOL isrow)
1645 {
1646 partialrec *blockdata;
1647
1648 blockdata = IF(isrow, lp->rowblocks, lp->colblocks);
1649 if(blockdata == NULL)
1650 return( 1 );
1651 else {
1652 if((blockdata->blocknow < 1) || (blockdata->blocknow > blockdata->blockcount))
1653 blockdata->blocknow = 1;
1654 return( blockdata->blockend[blockdata->blocknow-1] );
1655 }
1656 }
partial_blockEnd(lprec * lp,MYBOOL isrow)1657 STATIC int partial_blockEnd(lprec *lp, MYBOOL isrow)
1658 {
1659 partialrec *blockdata;
1660
1661 blockdata = IF(isrow, lp->rowblocks, lp->colblocks);
1662 if(blockdata == NULL)
1663 return( IF(isrow, lp->rows, lp->sum) );
1664 else {
1665 if((blockdata->blocknow < 1) || (blockdata->blocknow > blockdata->blockcount))
1666 blockdata->blocknow = 1;
1667 return( blockdata->blockend[blockdata->blocknow]-1 );
1668 }
1669 }
partial_blockNextPos(lprec * lp,int block,MYBOOL isrow)1670 STATIC int partial_blockNextPos(lprec *lp, int block, MYBOOL isrow)
1671 {
1672 partialrec *blockdata;
1673
1674 blockdata = IF(isrow, lp->rowblocks, lp->colblocks);
1675 #ifdef Paranoia
1676 if((blockdata == NULL) || (block <= 1) || (block > blockdata->blockcount)) {
1677 report(lp, SEVERE, "partial_blockNextPos: Invalid block %d specified.\n",
1678 block);
1679 return( -1 );
1680 }
1681 #endif
1682 block--;
1683 if(blockdata->blockpos[block] == blockdata->blockend[block+1])
1684 blockdata->blockpos[block] = blockdata->blockend[block];
1685 else
1686 blockdata->blockpos[block]++;
1687 return( blockdata->blockpos[block] );
1688 }
partial_blockStep(lprec * lp,MYBOOL isrow)1689 STATIC MYBOOL partial_blockStep(lprec *lp, MYBOOL isrow)
1690 {
1691 partialrec *blockdata;
1692
1693 blockdata = IF(isrow, lp->rowblocks, lp->colblocks);
1694 if(blockdata == NULL)
1695 return( FALSE );
1696 else if(blockdata->blocknow < blockdata->blockcount) {
1697 blockdata->blocknow++;
1698 return( TRUE);
1699 }
1700 else {
1701 blockdata->blocknow = 1;
1702 return( TRUE );
1703 }
1704 }
partial_isVarActive(lprec * lp,int varno,MYBOOL isrow)1705 STATIC MYBOOL partial_isVarActive(lprec *lp, int varno, MYBOOL isrow)
1706 {
1707 partialrec *blockdata;
1708
1709 blockdata = IF(isrow, lp->rowblocks, lp->colblocks);
1710 if(blockdata == NULL)
1711 return( TRUE );
1712 else {
1713 return( (MYBOOL) ((varno >= blockdata->blockend[blockdata->blocknow-1]) &&
1714 (varno < blockdata->blockend[blockdata->blocknow])) );
1715 }
1716 }
1717
1718
1719 /* Multiple pricing routines */
multi_create(lprec * lp,MYBOOL truncinf)1720 STATIC multirec *multi_create(lprec *lp, MYBOOL truncinf)
1721 {
1722 multirec *multi;
1723
1724 multi = (multirec *) calloc(1, sizeof(*multi));
1725 if(multi != NULL) {
1726 multi->active = 1;
1727 multi->lp = lp;
1728 multi->epszero = lp->epsprimal;
1729 multi->truncinf = truncinf;
1730 }
1731
1732 return(multi);
1733 }
multi_free(multirec ** multi)1734 STATIC void multi_free(multirec **multi)
1735 {
1736 if((multi == NULL) || (*multi == NULL))
1737 return;
1738 FREE((*multi)->items);
1739 FREE((*multi)->valueList);
1740 FREE((*multi)->indexSet);
1741 FREE((*multi)->freeList);
1742 FREE((*multi)->sortedList);
1743 FREE(*multi);
1744 }
multi_mustupdate(multirec * multi)1745 STATIC MYBOOL multi_mustupdate(multirec *multi)
1746 {
1747 return( (MYBOOL) ((multi != NULL) &&
1748 (multi->used < multi->limit)) );
1749 }
multi_resize(multirec * multi,int blocksize,int blockdiv,MYBOOL doVlist,MYBOOL doIset)1750 STATIC MYBOOL multi_resize(multirec *multi, int blocksize, int blockdiv, MYBOOL doVlist, MYBOOL doIset)
1751 {
1752 MYBOOL ok = TRUE;
1753
1754 if((blocksize > 1) && (blockdiv > 0)) {
1755 int oldsize = multi->size;
1756
1757 multi->size = blocksize;
1758 if(blockdiv > 1)
1759 multi->limit += (multi->size-oldsize) / blockdiv;
1760
1761 multi->items = (pricerec *) realloc(multi->items, (multi->size+1)*sizeof(*(multi->items)));
1762 multi->sortedList = (UNIONTYPE QSORTrec *) realloc(multi->sortedList, (multi->size+1)*sizeof(*(multi->sortedList)));
1763 ok = (multi->items != NULL) && (multi->sortedList != NULL) &&
1764 allocINT(multi->lp, &(multi->freeList), multi->size+1, AUTOMATIC);
1765 if(ok) {
1766 int i, n;
1767
1768 if(oldsize == 0)
1769 i = 0;
1770 else
1771 i = multi->freeList[0];
1772 multi->freeList[0] = i + (multi->size-oldsize);
1773 for(n = multi->size - 1, i++; i <= multi->freeList[0]; i++, n--)
1774 multi->freeList[i] = n;
1775 }
1776 if(doVlist)
1777 ok &= allocREAL(multi->lp, &(multi->valueList), multi->size+1, AUTOMATIC);
1778 if(doIset) {
1779 ok &= allocINT(multi->lp, &(multi->indexSet), multi->size+1, AUTOMATIC);
1780 if(ok && (oldsize == 0))
1781 multi->indexSet[0] = 0;
1782 }
1783 if(!ok)
1784 goto Undo;
1785
1786 }
1787 else {
1788 Undo:
1789 multi->size = 0;
1790 FREE(multi->items);
1791 FREE(multi->valueList);
1792 FREE(multi->indexSet);
1793 FREE(multi->freeList);
1794 FREE(multi->sortedList);
1795 }
1796 multi->active = 1;
1797
1798 return( ok );
1799 }
1800
multi_size(multirec * multi)1801 STATIC int multi_size(multirec *multi)
1802 {
1803 if(multi == NULL)
1804 return( 0 );
1805 else
1806 return( multi->size );
1807 }
1808
multi_used(multirec * multi)1809 STATIC int multi_used(multirec *multi)
1810 {
1811 if(multi == NULL)
1812 return( 0 );
1813 else
1814 return( multi->used );
1815 }
1816
multi_restart(multirec * multi)1817 STATIC int multi_restart(multirec *multi)
1818 {
1819 int i, n = multi->used;
1820
1821 multi->used = 0;
1822 multi->sorted = FALSE;
1823 multi->dirty = FALSE;
1824 if(multi->freeList != NULL) {
1825 for(i = 1; i <= multi->size; i++)
1826 multi->freeList[i] = multi->size - i;
1827 multi->freeList[0] = multi->size;
1828 }
1829 #if 0
1830 if(multi->indexSet != NULL)
1831 multi->indexSet[0] = 0;
1832 #endif
1833 return( n );
1834 }
1835
multi_valueInit(multirec * multi,REAL step_base,REAL obj_base)1836 STATIC void multi_valueInit(multirec *multi, REAL step_base, REAL obj_base)
1837 {
1838 multi->step_base = multi->step_last = step_base;
1839 multi->obj_base = multi->obj_last = obj_base;
1840 #ifdef Paranoia
1841 if(step_base > 0)
1842 report(multi->lp, SEVERE, "multi_valueInit: Positive constraint violation %g provided at iteration %6.0f\n",
1843 step_base, (double) get_total_iter(multi->lp));
1844 #endif
1845 }
1846
multi_valueList(multirec * multi)1847 STATIC REAL *multi_valueList(multirec *multi)
1848 {
1849 return(multi->valueList);
1850 }
1851
multi_indexSet(multirec * multi,MYBOOL regenerate)1852 STATIC int *multi_indexSet(multirec *multi, MYBOOL regenerate)
1853 {
1854 if(regenerate)
1855 multi_populateSet(multi, NULL, -1);
1856 return(multi->indexSet);
1857 }
1858
multi_getvar(multirec * multi,int item)1859 STATIC int multi_getvar(multirec *multi, int item)
1860 {
1861 #ifdef Paranoia
1862 if((item < 1) || (item >= multi->size))
1863 return(-1);
1864 #endif
1865 return( ((pricerec *) &(multi->sortedList[item].pvoidreal.ptr))->varno );
1866 }
1867
multi_recompute(multirec * multi,int index,MYBOOL isphase2,MYBOOL fullupdate)1868 STATIC MYBOOL multi_recompute(multirec *multi, int index, MYBOOL isphase2, MYBOOL fullupdate)
1869 {
1870 int i, n;
1871 REAL lB, uB, Alpha, this_theta, prev_theta;
1872 lprec *lp = multi->lp;
1873 pricerec *thisprice;
1874
1875 /* Define target update window */
1876 if(multi->dirty) {
1877 index = 0;
1878 n = multi->used - 1;
1879 }
1880 else if(fullupdate)
1881 n = multi->used - 1;
1882 else
1883 n = index;
1884
1885 /* Initialize accumulators from the specified update index */
1886 if(index == 0) {
1887 multi->maxpivot = 0;
1888 multi->maxbound = 0;
1889 multi->step_last = multi->step_base;
1890 multi->obj_last = multi->obj_base;
1891 thisprice = NULL;
1892 this_theta = 0;
1893 }
1894 else {
1895 multi->obj_last = multi->valueList[index-1];
1896 multi->step_last = multi->sortedList[index-1].pvoidreal.realval;
1897 thisprice = (pricerec *) (multi->sortedList[index-1].pvoidreal.ptr);
1898 this_theta = thisprice->theta;
1899 }
1900
1901 /* Update step lengths and objective values */
1902 while((index <= n) && (multi->step_last < multi->epszero)) {
1903
1904 /* Update parameters for this loop */
1905 prev_theta = this_theta;
1906 thisprice = (pricerec *) (multi->sortedList[index].pvoidreal.ptr);
1907 this_theta = thisprice->theta;
1908 Alpha = fabs(thisprice->pivot);
1909 uB = lp->upbo[thisprice->varno];
1910 lB = 0;
1911 SETMAX(multi->maxpivot, Alpha);
1912 SETMAX(multi->maxbound, uB);
1913
1914 /* Do the value updates */
1915 if(isphase2) {
1916 multi->obj_last += (this_theta - prev_theta) * multi->step_last; /* Sign-readjusted from coldual()/Maros */
1917 if(uB >= lp->infinite)
1918 multi->step_last = lp->infinite;
1919 else
1920 multi->step_last += Alpha*(uB-lB);
1921 }
1922 else {
1923 multi->obj_last += (this_theta - prev_theta) * multi->step_last; /* Sign-readjusted from coldual()/Maros */
1924 multi->step_last += Alpha;
1925 }
1926
1927 /* Store updated values at the indexed locations */
1928 multi->sortedList[index].pvoidreal.realval = multi->step_last;
1929 multi->valueList[index] = multi->obj_last;
1930 #ifdef Paranoia
1931 if(lp->spx_trace &&
1932 (multi->step_last > lp->infinite))
1933 report(lp, SEVERE, "multi_recompute: A very large step-size %g was generated at iteration %6.0f\n",
1934 multi->step_last, (double) get_total_iter(lp));
1935 #endif
1936 index++;
1937 }
1938
1939 /* Discard candidates entered earlier that now make the OF worsen, and
1940 make sure that the released positions are added to the free list. */
1941 n = index;
1942 while(n < multi->used) {
1943 i = ++multi->freeList[0];
1944 multi->freeList[i] = ((pricerec *) multi->sortedList[n].pvoidreal.ptr) - multi->items;
1945 n++;
1946 }
1947 multi->used = index;
1948 if(multi->sorted && (index == 1))
1949 multi->sorted = FALSE;
1950 multi->dirty = FALSE;
1951
1952 /* Return TRUE if the step is now positive */
1953 return( (MYBOOL) (multi->step_last >= multi->epszero) );
1954 }
1955
multi_truncatingvar(multirec * multi,int varnr)1956 STATIC MYBOOL multi_truncatingvar(multirec *multi, int varnr)
1957 {
1958 return( multi->truncinf && is_infinite(multi->lp, multi->lp->upbo[varnr]) );
1959 }
1960
multi_removevar(multirec * multi,int varnr)1961 STATIC MYBOOL multi_removevar(multirec *multi, int varnr)
1962 {
1963 int i = 1;
1964 int *coltarget = multi->indexSet;
1965
1966 if(coltarget == NULL)
1967 return( FALSE );
1968
1969 while((i <= multi->used) && (coltarget[i] != varnr))
1970 i++;
1971 if(i > multi->used)
1972 return( FALSE );
1973
1974 for(; i < multi->used; i++)
1975 coltarget[i] = coltarget[i+1];
1976 coltarget[0]--;
1977 multi->used--;
1978 multi->dirty = TRUE;
1979 return( TRUE );
1980 }
1981
multi_enteringvar(multirec * multi,pricerec * current,int priority)1982 STATIC int multi_enteringvar(multirec *multi, pricerec *current, int priority)
1983 {
1984 lprec *lp = multi->lp;
1985 int i, bestindex, colnr;
1986 REAL bound, score, bestscore = -lp->infinite;
1987 REAL b1, b2, b3;
1988 pricerec *candidate, *bestcand;
1989
1990 i = 0;
1991
1992 /* Check that we have a candidate */
1993 multi->active = bestindex = 0;
1994 if((multi == NULL) || (multi->used == 0))
1995 return( bestindex );
1996
1997 /* Check for pruning possibility of the B&B tree */
1998 if(multi->objcheck && (lp->solutioncount > 0) &&
1999 bb_better(lp, OF_WORKING | OF_PROJECTED, OF_TEST_WE)) {
2000 lp->spx_status = FATHOMED;
2001 return( bestindex );
2002 }
2003
2004 /* Check the trivial case */
2005 if(multi->used == 1) {
2006 bestcand = (pricerec *) (multi->sortedList[bestindex].pvoidreal.ptr);
2007 goto Finish;
2008 }
2009
2010 /* Set priority weights */
2011 Redo:
2012 switch(priority) {
2013 case 0: b1 = 0.0, b2 = 0.0, b3 = 1.0; /* Only OF */
2014 bestindex = multi->used - 2; break;
2015 case 1: b1 = 0.2, b2 = 0.3, b3 = 0.5; break; /* Emphasize OF */
2016 case 2: b1 = 0.3, b2 = 0.5, b3 = 0.2; break; /* Emphasize bound */
2017 case 3: b1 = 0.6, b2 = 0.2, b3 = 0.2; break; /* Emphasize pivot */
2018 case 4: b1 = 1.0, b2 = 0.0, b3 = 0.0; break; /* Only pivot */
2019 default: b1 = 0.4, b2 = 0.2, b3 = 0.4; /* Balanced default */
2020 }
2021 bestcand = (pricerec *) (multi->sortedList[bestindex].pvoidreal.ptr);
2022
2023 /* Loop over all candidates to get the best entering candidate;
2024 start at the end to try to maximize the chain length */
2025 for(i = multi->used - 1; i >= 0; i--) {
2026 candidate = (pricerec *) (multi->sortedList[i].pvoidreal.ptr);
2027 colnr = candidate->varno;
2028 bound = lp->upbo[colnr];
2029 score = fabs(candidate->pivot) / multi->maxpivot;
2030 score = pow(1.0 + score , b1) *
2031 pow(1.0 + log(bound / multi->maxbound + 1), b2) *
2032 pow(1.0 + (REAL) i / multi->used , b3);
2033 if(score > bestscore) {
2034 bestscore = score;
2035 bestindex = i;
2036 bestcand = candidate;
2037 }
2038 }
2039
2040 /* Do pivot protection */
2041 if((priority < 4) && (fabs(bestcand->pivot) < lp->epssolution)) {
2042 bestindex = 0;
2043 priority++;
2044 goto Redo;
2045 }
2046
2047 Finish:
2048 /* Make sure we shrink the list and update */
2049 multi->active = colnr = bestcand->varno;
2050 if(bestindex < multi->used - 1) {
2051 #if 0
2052 /* if(lp->upbo[colnr] >= lp->infinite) */
2053 QS_swap(multi->sortedList, bestindex, multi->used-1);
2054 multi_recompute(multi, bestindex, (bestcand->isdual == AUTOMATIC), TRUE);
2055 #else
2056 multi->used = i + 1;
2057 #endif
2058 }
2059 multi_populateSet(multi, NULL, multi->active);
2060
2061 /* Compute the entering theta and update parameters */
2062 score = (multi->used == 1 ? multi->step_base : multi->sortedList[multi->used-2].pvoidreal.realval);
2063 score /= bestcand->pivot;
2064 score = my_chsign(!lp->is_lower[multi->active], score);
2065
2066 if(lp->spx_trace &&
2067 (fabs(score) > 1/lp->epsprimal))
2068 report(lp, IMPORTANT, "multi_enteringvar: A very large Theta %g was generated (pivot %g)\n",
2069 score, bestcand->pivot);
2070 multi->step_base = score;
2071 if(current != NULL)
2072 *current = *bestcand;
2073
2074 return( multi->active );
2075 }
2076
multi_enteringtheta(multirec * multi)2077 STATIC REAL multi_enteringtheta(multirec *multi)
2078 {
2079 return( multi->step_base );
2080 }
2081
multi_populateSet(multirec * multi,int ** list,int excludenr)2082 STATIC int multi_populateSet(multirec *multi, int **list, int excludenr)
2083 {
2084 int n = 0;
2085 if(list == NULL)
2086 list = &(multi->indexSet);
2087 if((multi->used > 0) &&
2088 ((*list != NULL) || allocINT(multi->lp, list, multi->size+1, FALSE))) {
2089 int i, colnr;
2090
2091 for(i = 0; i < multi->used; i++) {
2092 colnr = ((pricerec *) (multi->sortedList[i].pvoidreal.ptr))->varno;
2093 if((colnr != excludenr) &&
2094 /* Prevent an unbounded variable from "bound-flip"; this could
2095 actually indicate that we should let the entering variable be
2096 bound-swapped (in the case that it is bounded), but we
2097 disregard this possibility here, since it brings with it
2098 issues of pivot size, etc. */
2099 ((excludenr > 0) && (multi->lp->upbo[colnr] < multi->lp->infinite))) {
2100 n++;
2101 (*list)[n] = colnr;
2102 }
2103 }
2104 (*list)[0] = n;
2105 }
2106 return( n );
2107 }
2108
2109