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(&current, &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(&current, &candidate, NULL))
1108           break;
1109         if((Hpass == 2) && (current.varno == candidate.varno))
1110           Htheta = f;
1111       }
1112       else
1113         if(findSubstitutionVar(&current, &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(&current, &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(&current, &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