1 
2 #include <string.h>
3 #include "commonlib.h"
4 #include "lp_lib.h"
5 #include "lp_report.h"
6 #include "lp_scale.h"
7 
8 #ifdef FORTIFY
9 # include "lp_fortify.h"
10 #endif
11 
12 
13 /*
14     Scaling routines for lp_solve v5.0+
15    ----------------------------------------------------------------------------------
16     Author:        Kjell Eikland
17     Contact:       kjell.eikland@broadpark.no
18     License terms: LGPL.
19 
20     Requires:      lp_lib.h, lp_scale.h
21 
22     Release notes:
23     v5.0.0  1 January 2004      Significantly expanded and repackaged scaling
24                                 routines.
25     v5.0.1  20 February 2004    Modified rounding behaviour in several areas.
26     v5.1.0  20 July 2004        Reworked with flexible matrix storage model.
27     v5.2.0  20 February 2005    Converted to matrix storage model without the OF.
28 
29    ----------------------------------------------------------------------------------
30 */
31 
32 /* First define scaling and unscaling primitives */
33 
scaled_value(lprec * lp,REAL value,int index)34 REAL scaled_value(lprec *lp, REAL value, int index)
35 {
36   if(fabs(value) < lp->infinite) {
37     if(lp->scaling_used) {
38       if(index > lp->rows)
39         value /= lp->scalars[index];
40       else
41         value *= lp->scalars[index];
42     }
43   }
44   else
45     value = my_sign(value)*lp->infinite;
46   return(value);
47 }
48 
unscaled_value(lprec * lp,REAL value,int index)49 REAL unscaled_value(lprec *lp, REAL value, int index)
50 {
51   if(fabs(value) < lp->infinite) {
52     if(lp->scaling_used) {
53       if(index > lp->rows)
54         value *= lp->scalars[index];
55       else
56         value /= lp->scalars[index];
57     }
58   }
59   else
60     value = my_sign(value)*lp->infinite;
61   return(value);
62 }
63 
scaled_mat(lprec * lp,REAL value,int rownr,int colnr)64 STATIC REAL scaled_mat(lprec *lp, REAL value, int rownr, int colnr)
65 {
66   if(lp->scaling_used)
67     value *= lp->scalars[rownr] * lp->scalars[lp->rows + colnr];
68   return( value );
69 }
70 
unscaled_mat(lprec * lp,REAL value,int rownr,int colnr)71 STATIC REAL unscaled_mat(lprec *lp, REAL value, int rownr, int colnr)
72 {
73   if(lp->scaling_used)
74     value /= lp->scalars[rownr] * lp->scalars[lp->rows + colnr];
75   return( value );
76 }
77 
78 /* Compute the scale factor by the formulae:
79       FALSE: SUM (log |Aij|) ^ 2
80       TRUE:  SUM (log |Aij| - RowScale[i] - ColScale[j]) ^ 2 */
CurtisReidMeasure(lprec * lp,MYBOOL _Advanced,REAL * FRowScale,REAL * FColScale)81 REAL CurtisReidMeasure(lprec *lp, MYBOOL _Advanced, REAL *FRowScale, REAL *FColScale)
82 {
83   int      i, nz;
84   REAL     absvalue, logvalue;
85   register REAL result;
86   MATrec   *mat = lp->matA;
87   REAL     *value;
88   int      *rownr, *colnr;
89 
90   /* Do OF part */
91   result = 0;
92   for(i = 1; i <= lp->columns; i++) {
93     absvalue = fabs(lp->orig_obj[i]);
94     if(absvalue > 0) {
95       logvalue = log(absvalue);
96       if(_Advanced)
97         logvalue -= FRowScale[0] + FColScale[i];
98       result += logvalue*logvalue;
99     }
100   }
101 
102   /* Do constraint matrix part */
103   mat_validate(mat);
104   value = &(COL_MAT_VALUE(0));
105   rownr = &(COL_MAT_ROWNR(0));
106   colnr = &(COL_MAT_COLNR(0));
107   nz = get_nonzeros(lp);
108   for(i = 0; i < nz;
109       i++, value += matValueStep, rownr += matRowColStep, colnr += matRowColStep) {
110     absvalue = fabs(*value);
111     if(absvalue > 0) {
112       logvalue = log(absvalue);
113       if(_Advanced)
114         logvalue -= FRowScale[*rownr] + FColScale[*colnr];
115       result += logvalue*logvalue;
116     }
117   }
118   return( result );
119 }
120 
121 /* Implementation of the Curtis-Reid scaling based on the paper
122    "On the Automatic Scaling of Matrices for Gaussian
123     Elimination," Journal of the Institute of Mathematics and
124     Its Applications (1972) 10, 118-124.
125 
126     Solve the system | M   E | (r)   (sigma)
127                      |       | ( ) = (     )
128                      | E^T N | (c)   ( tau )
129 
130     by the conjugate gradient method (clever recurrences).
131 
132     E is the matrix A with all elements = 1
133 
134     M is diagonal matrix of row    counts (RowCount)
135     N is diagonal matrix of column counts (ColCount)
136 
137     sigma is the vector of row    logarithm sums (RowSum)
138     tau   is the vector of column logarithm sums (ColSum)
139 
140     r, c are resulting row and column scalings (RowScale, ColScale) */
141 
CurtisReidScales(lprec * lp,MYBOOL _Advanced,REAL * FRowScale,REAL * FColScale)142 int CurtisReidScales(lprec *lp, MYBOOL _Advanced, REAL *FRowScale, REAL *FColScale)
143 {
144   int    i, row, col, ent, nz;
145   REAL   *RowScalem2, *ColScalem2,
146          *RowSum, *ColSum,
147          *residual_even, *residual_odd;
148   REAL   sk,   qk,     ek,
149          skm1, qkm1,   ekm1,
150          qkm2, qkqkm1, ekm2, ekekm1,
151          absvalue, logvalue,
152          StopTolerance;
153   int    *RowCount, *ColCount, colMax;
154   int    Result;
155   MATrec *mat = lp->matA;
156   REAL   *value;
157   int    *rownr, *colnr;
158 
159   if(CurtisReidMeasure(lp, _Advanced, FRowScale, FColScale)<0.1*get_nonzeros(lp))
160   return(0);
161 
162   /* Allocate temporary memory and find RowSum and ColSum measures */
163   nz = get_nonzeros(lp);
164   colMax = lp->columns;
165 
166   allocREAL(lp, &RowSum, lp->rows+1, TRUE);
167   allocINT(lp,  &RowCount, lp->rows+1, TRUE);
168   allocREAL(lp, &residual_odd, lp->rows+1, TRUE);
169 
170   allocREAL(lp, &ColSum, colMax+1, TRUE);
171   allocINT(lp,  &ColCount, colMax+1, TRUE);
172   allocREAL(lp, &residual_even, colMax+1, TRUE);
173 
174   allocREAL(lp, &RowScalem2, lp->rows+1, FALSE);
175   allocREAL(lp, &ColScalem2, colMax+1, FALSE);
176 
177   /* Set origin for row scaling */
178   for(i = 1; i <= colMax; i++) {
179     absvalue=fabs(lp->orig_obj[i]);
180     if(absvalue>0) {
181       logvalue = log(absvalue);
182       ColSum[i] += logvalue;
183       RowSum[0] += logvalue;
184       ColCount[i]++;
185       RowCount[0]++;
186     }
187   }
188 
189   value = &(COL_MAT_VALUE(0));
190   rownr = &(COL_MAT_ROWNR(0));
191   colnr = &(COL_MAT_COLNR(0));
192   for(i = 0; i < nz;
193       i++, value += matValueStep, rownr += matRowColStep, colnr += matRowColStep) {
194     absvalue=fabs(*value);
195     if(absvalue>0) {
196       logvalue = log(absvalue);
197       ColSum[*colnr] += logvalue;
198       RowSum[*rownr] += logvalue;
199       ColCount[*colnr]++;
200       RowCount[*rownr]++;
201     }
202   }
203 
204   /* Make sure we dont't have division by zero errors */
205   for(row = 0; row <= lp->rows; row++)
206     if(RowCount[row] == 0)
207       RowCount[row] = 1;
208   for(col = 1; col <= colMax; col++)
209     if(ColCount[col] == 0)
210       ColCount[col] = 1;
211 
212   /* Initialize to RowScale = RowCount-1 RowSum
213                    ColScale = 0.0
214                    residual = ColSum - ET RowCount-1 RowSum */
215 
216   StopTolerance= MAX(lp->scalelimit-floor(lp->scalelimit), DEF_SCALINGEPS);
217   StopTolerance *= (REAL) nz;
218   for(row = 0; row <= lp->rows; row++) {
219     FRowScale[row] = RowSum[row] / (REAL) RowCount[row];
220     RowScalem2[row] = FRowScale[row];
221   }
222 
223   /* Compute initial residual */
224   for(col = 1; col <= colMax; col++) {
225     FColScale[col] = 0;
226     ColScalem2[col] = 0;
227     residual_even[col] = ColSum[col];
228 
229     if(lp->orig_obj[col] != 0)
230       residual_even[col] -= RowSum[0] / (REAL) RowCount[0];
231 
232     i = mat->col_end[col-1];
233     rownr = &(COL_MAT_ROWNR(i));
234     ent = mat->col_end[col];
235     for(; i < ent;
236         i++, rownr += matRowColStep) {
237       residual_even[col] -= RowSum[*rownr] / (REAL) RowCount[*rownr];
238     }
239   }
240 
241   /* Compute sk */
242   sk = 0;
243   skm1 = 0;
244   for(col = 1; col <= colMax; col++)
245     sk += (residual_even[col]*residual_even[col]) / (REAL) ColCount[col];
246 
247   Result = 0;
248   qk=1; qkm1=0; qkm2=0;
249   ek=0; ekm1=0; ekm2=0;
250 
251   while(sk>StopTolerance) {
252   /* Given the values of residual and sk, construct
253      ColScale (when pass is even)
254      RowScale (when pass is odd)  */
255 
256     qkqkm1 = qk * qkm1;
257     ekekm1 = ek * ekm1;
258     if((Result % 2) == 0) { /* pass is even; construct RowScale[pass+1] */
259       if(Result != 0) {
260         for(row = 0; row <= lp->rows; row++)
261           RowScalem2[row] = FRowScale[row];
262         if(qkqkm1 != 0) {
263           for(row = 0; row <= lp->rows; row++)
264             FRowScale[row]*=(1 + ekekm1 / qkqkm1);
265           for(row = 0; row<=lp->rows; row++)
266             FRowScale[row]+=(residual_odd[row] / (qkqkm1 * (REAL) RowCount[row]) -
267                              RowScalem2[row] * ekekm1 / qkqkm1);
268         }
269       }
270     }
271     else { /* pass is odd; construct ColScale[pass+1] */
272       for(col = 1; col <= colMax; col++)
273         ColScalem2[col] = FColScale[col];
274       if(qkqkm1 != 0) {
275         for(col = 1; col <= colMax; col++)
276           FColScale[col] *= (1 + ekekm1 / qkqkm1);
277         for(col = 1; col <= colMax; col++)
278           FColScale[col] += (residual_even[col] / ((REAL) ColCount[col] * qkqkm1) -
279                              ColScalem2[col] * ekekm1 / qkqkm1);
280       }
281     }
282 
283     /* update residual and sk (pass + 1) */
284     if((Result % 2) == 0) { /* even */
285        /* residual */
286       for(row = 0; row <= lp->rows; row++)
287         residual_odd[row] *= ek;
288 
289       for(i = 1; i <= colMax; i++)
290         if(lp->orig_obj[i] != 0)
291           residual_odd[0] += (residual_even[i] / (REAL) ColCount[i]);
292 
293       rownr = &(COL_MAT_ROWNR(0));
294       colnr = &(COL_MAT_COLNR(0));
295       for(i = 0; i < nz;
296           i++, rownr += matRowColStep, colnr += matRowColStep) {
297         residual_odd[*rownr] += (residual_even[*colnr] / (REAL) ColCount[*colnr]);
298       }
299       for(row = 0; row <= lp->rows; row++)
300         residual_odd[row] *= (-1 / qk);
301 
302       /* sk */
303       skm1 = sk;
304       sk = 0;
305       for(row = 0; row <= lp->rows; row++)
306         sk += (residual_odd[row]*residual_odd[row]) / (REAL) RowCount[row];
307     }
308     else { /* odd */
309       /* residual */
310       for(col = 1; col <= colMax; col++)
311         residual_even[col] *= ek;
312 
313       for(i = 1; i <= colMax; i++)
314         if(lp->orig_obj[i] != 0)
315           residual_even[i] += (residual_odd[0] / (REAL) RowCount[0]);
316 
317       rownr = &(COL_MAT_ROWNR(0));
318       colnr = &(COL_MAT_COLNR(0));
319       for(i = 0; i < nz;
320           i++, rownr += matRowColStep, colnr += matRowColStep) {
321         residual_even[*colnr] += (residual_odd[*rownr] / (REAL) RowCount[*rownr]);
322       }
323       for(col = 1; col <= colMax; col++)
324         residual_even[col] *= (-1 / qk);
325 
326       /* sk */
327       skm1 = sk;
328       sk = 0;
329       for(col = 1; col <= colMax; col++)
330         sk += (residual_even[col]*residual_even[col]) / (REAL) ColCount[col];
331     }
332 
333     /* Compute ek and qk */
334     ekm2=ekm1;
335     ekm1=ek;
336     ek=qk * sk / skm1;
337 
338     qkm2=qkm1;
339     qkm1=qk;
340     qk=1-ek;
341 
342     Result++;
343   }
344 
345   /* Synchronize the RowScale and ColScale vectors */
346   ekekm1 = ek * ekm1;
347   if(qkm1 != 0) {
348   if((Result % 2) == 0) { /* pass is even, compute RowScale */
349     for(row = 0; row<=lp->rows; row++)
350       FRowScale[row]*=(1.0 + ekekm1 / qkm1);
351     for(row = 0; row<=lp->rows; row++)
352       FRowScale[row]+=(residual_odd[row] / (qkm1 * (REAL) RowCount[row]) -
353                       RowScalem2[row] * ekekm1 / qkm1);
354     }
355   }
356   else { /* pass is odd, compute ColScale */
357     for(col=1; col<=colMax; col++)
358       FColScale[col]*=(1 + ekekm1 / qkm1);
359     for(col=1; col<=colMax; col++)
360       FColScale[col]+=(residual_even[col] / ((REAL) ColCount[col] * qkm1) -
361                        ColScalem2[col] * ekekm1 / qkm1);
362   }
363 
364   /* Do validation, if indicated */
365   if(FALSE && mat_validate(mat)){
366     double check, error;
367 
368     /* CHECK: M RowScale + E ColScale = RowSum */
369     error = 0;
370     for(row = 0; row <= lp->rows; row++) {
371       check = (REAL) RowCount[row] * FRowScale[row];
372       if(row == 0) {
373         for(i = 1; i <= colMax; i++) {
374           if(lp->orig_obj[i] != 0)
375             check += FColScale[i];
376         }
377       }
378       else {
379         i = mat->row_end[row-1];
380         ent = mat->row_end[row];
381         for(; i < ent; i++) {
382           col = ROW_MAT_COLNR(i);
383           check += FColScale[col];
384         }
385       }
386       check -= RowSum[row];
387       error += check*check;
388     }
389 
390     /* CHECK: E^T RowScale + N ColScale = ColSum */
391     error = 0;
392     for(col = 1; col <= colMax; col++) {
393       check = (REAL) ColCount[col] * FColScale[col];
394 
395       if(lp->orig_obj[col] != 0)
396         check += FRowScale[0];
397 
398       i = mat->col_end[col-1];
399       ent = mat->col_end[col];
400       rownr = &(COL_MAT_ROWNR(i));
401       for(; i < ent;
402           i++, rownr += matRowColStep) {
403         check += FRowScale[*rownr];
404       }
405       check -= ColSum[col];
406       error += check*check;
407     }
408   }
409 
410   /* Convert to scaling factors (rounding to nearest power
411      of 2 can optionally be done as a separate step later) */
412   for(col = 1; col <= colMax; col++) {
413     absvalue = exp(-FColScale[col]);
414     if(absvalue < MIN_SCALAR) absvalue = MIN_SCALAR;
415     if(absvalue > MAX_SCALAR) absvalue = MAX_SCALAR;
416     if(!is_int(lp,col) || is_integerscaling(lp))
417         FColScale[col] = absvalue;
418     else
419         FColScale[col] = 1;
420   }
421   for(row = 0; row <= lp->rows; row++) {
422     absvalue = exp(-FRowScale[row]);
423     if(absvalue < MIN_SCALAR) absvalue = MIN_SCALAR;
424     if(absvalue > MAX_SCALAR) absvalue = MAX_SCALAR;
425     FRowScale[row] = absvalue;
426   }
427 
428  /* free temporary memory */
429   FREE(RowSum);
430   FREE(ColSum);
431   FREE(RowCount);
432   FREE(ColCount);
433   FREE(residual_even);
434   FREE(residual_odd);
435   FREE(RowScalem2);
436   FREE(ColScalem2);
437 
438   return(Result);
439 
440 }
441 
scaleCR(lprec * lp,REAL * scaledelta)442 STATIC MYBOOL scaleCR(lprec *lp, REAL *scaledelta)
443 {
444   REAL *scalechange = NULL;
445   int  Result;
446 
447   if(!lp->scaling_used) {
448     allocREAL(lp, &lp->scalars, lp->sum_alloc + 1, FALSE);
449     for(Result = 0; Result <= lp->sum; Result++)
450       lp->scalars[Result] = 1;
451     lp->scaling_used = TRUE;
452   }
453 
454   if(scaledelta == NULL)
455     allocREAL(lp, &scalechange, lp->sum + 1, FALSE);
456   else
457     scalechange = scaledelta;
458 
459   Result=CurtisReidScales(lp, FALSE, scalechange, &scalechange[lp->rows]);
460   if(Result>0) {
461 
462     /* Do the scaling*/
463     if(scale_updaterows(lp, scalechange, TRUE) ||
464        scale_updatecolumns(lp, &scalechange[lp->rows], TRUE))
465       lp->scalemode |= SCALE_CURTISREID;
466 
467     set_action(&lp->spx_action, ACTION_REBASE | ACTION_REINVERT | ACTION_RECOMPUTE);
468   }
469 
470   if(scaledelta == NULL)
471     FREE(scalechange);
472 
473   return((MYBOOL) (Result > 0));
474 }
475 
transform_for_scale(lprec * lp,REAL * value)476 STATIC MYBOOL transform_for_scale(lprec *lp, REAL *value)
477 {
478   MYBOOL Accept = TRUE;
479   *value = fabs(*value);
480 #ifdef Paranoia
481   if(*value < lp->epsmachine) {
482     Accept = FALSE;
483     report(lp, SEVERE, "transform_for_scale: A zero-valued entry was passed\n");
484   }
485   else
486 #endif
487   if(is_scalemode(lp, SCALE_LOGARITHMIC))
488     *value = log(*value);
489   else if(is_scalemode(lp, SCALE_QUADRATIC))
490     (*value) *= (*value);
491   return( Accept );
492 }
493 
accumulate_for_scale(lprec * lp,REAL * min,REAL * max,REAL value)494 STATIC void accumulate_for_scale(lprec *lp, REAL *min, REAL *max, REAL value)
495 {
496   if(transform_for_scale(lp, &value)) {
497     if(is_scaletype(lp, SCALE_MEAN)) {
498       *max += value;
499       *min += 1;
500     }
501     else {
502       SETMAX(*max, value);
503       SETMIN(*min, value);
504     }
505    }
506 }
507 
minmax_to_scale(lprec * lp,REAL min,REAL max,int itemcount)508 STATIC REAL minmax_to_scale(lprec *lp, REAL min, REAL max, int itemcount)
509 {
510   REAL scale;
511 
512   /* Initialize according to transformation / weighting model */
513   if(is_scalemode(lp, SCALE_LOGARITHMIC))
514     scale = 0;
515   else
516     scale = 1;
517   if(itemcount <= 0)
518     return(scale);
519 
520   /* Compute base scalar according to chosen scaling type */
521   if(is_scaletype(lp, SCALE_MEAN)) {
522     if(min > 0)
523       scale = max / min;
524   }
525   else if(is_scaletype(lp, SCALE_RANGE))
526     scale = (max + min) / 2;
527   else if(is_scaletype(lp, SCALE_GEOMETRIC))
528     scale = sqrt(min*max);
529   else if(is_scaletype(lp, SCALE_EXTREME))
530     scale = max;
531 
532   /* Compute final scalar according to transformation / weighting model */
533   if(is_scalemode(lp, SCALE_LOGARITHMIC))
534     scale = exp(-scale);
535   else if(is_scalemode(lp, SCALE_QUADRATIC)) {
536     if(scale == 0)
537       scale = 1;
538     else
539       scale = 1 / sqrt(scale);
540   }
541   else {
542     if(scale == 0)
543       scale = 1;
544     else
545       scale = 1 / scale;
546   }
547 
548   /* Make sure we are within acceptable scaling ranges */
549   SETMAX(scale, MIN_SCALAR);
550   SETMIN(scale, MAX_SCALAR);
551 
552   return(scale);
553 }
554 
roundPower2(REAL scale)555 STATIC REAL roundPower2(REAL scale)
556 /* Purpose is to round a number to it nearest power of 2; in a system
557    with binary number representation, this avoids rounding errors when
558    scale is used to normalize another value */
559 {
560   long int power2;
561   MYBOOL   isSmall = FALSE;
562 
563   if(scale == 1)
564     return( scale );
565 
566   /* Obtain the fractional power of 2 */
567   if(scale < 2) {
568     scale = 2 / scale;
569     isSmall = TRUE;
570   }
571   else
572     scale /= 2;
573   scale = log(scale)/log(2.0);
574 
575   /* Find the desired nearest power of two and compute the associated scalar */
576   power2 = (long) ceil(scale-0.5);
577   scale = 1 << power2;
578   if(isSmall)
579     scale = 1.0 / scale;
580 
581   return( scale );
582 
583 }
584 
scale_updatecolumns(lprec * lp,REAL * scalechange,MYBOOL updateonly)585 STATIC MYBOOL scale_updatecolumns(lprec *lp, REAL *scalechange, MYBOOL updateonly)
586 {
587   int i, j;
588 
589   /* Verify that the scale change is significant (different from the unit) */
590   for(i = lp->columns; i > 0; i--)
591     if(fabs(scalechange[i]-1) > lp->epsprimal)
592       break;
593   if(i <= 0)
594     return( FALSE );
595 
596  /* Update the pre-existing column scalar */
597   if(updateonly)
598     for(i = 1, j = lp->rows + 1; j <= lp->sum; i++, j++)
599       lp->scalars[j] *= scalechange[i];
600   else
601     for(i = 1, j = lp->rows + 1; j <= lp->sum; i++, j++)
602       lp->scalars[j] = scalechange[i];
603 
604   return( TRUE );
605 }
606 
scale_updaterows(lprec * lp,REAL * scalechange,MYBOOL updateonly)607 STATIC MYBOOL scale_updaterows(lprec *lp, REAL *scalechange, MYBOOL updateonly)
608 {
609   int i;
610 
611   /* Verify that the scale change is significant (different from the unit) */
612   for(i = lp->rows; i >= 0; i--) {
613     if(fabs(scalechange[i]-1) > lp->epsprimal)
614       break;
615   }
616   if(i < 0)
617     return( FALSE );
618 
619  /* Update the pre-existing row scalar */
620   if(updateonly)
621     for(i = 0; i <= lp->rows; i++)
622       lp->scalars[i] *= scalechange[i];
623   else
624     for(i = 0; i <= lp->rows; i++)
625       lp->scalars[i] = scalechange[i];
626 
627   return( TRUE );
628 }
629 
scale_columns(lprec * lp,REAL * scaledelta)630 STATIC MYBOOL scale_columns(lprec *lp, REAL *scaledelta)
631 {
632   int     i,j, colMax, nz;
633   REAL    *scalechange;
634   REAL    *value;
635   int     *colnr;
636   MATrec  *mat = lp->matA;
637 
638   /* Check that columns are in fact targeted */
639   if((lp->scalemode & SCALE_ROWSONLY) != 0)
640     return( TRUE );
641 
642   if(scaledelta == NULL)
643     scalechange = &lp->scalars[lp->rows];
644   else
645     scalechange = &scaledelta[lp->rows];
646 
647   colMax = lp->columns;
648 
649   /* Scale matrix entries (including any Lagrangean constraints) */
650   for(i = 1; i <= lp->columns; i++) {
651     lp->orig_obj[i] *= scalechange[i];
652   }
653 
654   mat_validate(lp->matA);
655   nz = get_nonzeros(lp);
656   value = &(COL_MAT_VALUE(0));
657   colnr = &(COL_MAT_COLNR(0));
658   for(i = 0; i < nz;
659       i++, value += matValueStep, colnr += matRowColStep) {
660     (*value) *= scalechange[*colnr];
661   }
662 
663   /* Scale variable bounds as well */
664   for(i = 1, j = lp->rows + 1; j <= lp->sum; i++, j++) {
665     if(lp->orig_lowbo[j] > -lp->infinite)
666       lp->orig_lowbo[j] /= scalechange[i];
667     if(lp->orig_upbo[j] < lp->infinite)
668       lp->orig_upbo[j] /= scalechange[i];
669     if(lp->sc_lobound[i] != 0)
670       lp->sc_lobound[i] /= scalechange[i];
671   }
672 
673   lp->columns_scaled = TRUE;
674   set_action(&lp->spx_action, ACTION_REBASE | ACTION_REINVERT | ACTION_RECOMPUTE);
675 
676   return( TRUE );
677 }
678 
scale_rows(lprec * lp,REAL * scaledelta)679 STATIC MYBOOL scale_rows(lprec *lp, REAL *scaledelta)
680 {
681   int     i, j, nz, colMax;
682   REAL    *scalechange;
683   REAL    *value;
684   int     *rownr;
685   MATrec  *mat = lp->matA;
686 
687 
688   /* Check that rows are in fact targeted */
689   if((lp->scalemode & SCALE_COLSONLY) != 0)
690     return( TRUE );
691 
692   if(scaledelta == NULL)
693     scalechange = lp->scalars;
694   else
695     scalechange = scaledelta;
696 
697   colMax = lp->columns;
698 
699   /* First row-scale the matrix (including the objective function) */
700   for(i = 1; i <= colMax; i++) {
701     lp->orig_obj[i] *= scalechange[0];
702   }
703 
704   nz = get_nonzeros(lp);
705   value = &(COL_MAT_VALUE(0));
706   rownr = &(COL_MAT_ROWNR(0));
707   for(i = 0; i < nz;
708       i++, value += matValueStep, rownr += matRowColStep) {
709     (*value) *= scalechange[*rownr];
710   }
711 
712   /* ...and scale the rhs and the row bounds (RANGES in MPS!!) */
713   for(i = 0; i <= lp->rows; i++) {
714     if(fabs(lp->orig_rhs[i]) < lp->infinite)
715       lp->orig_rhs[i] *= scalechange[i];
716 
717     j = lp->presolve_undo->var_to_orig[i];
718     if(j != 0)
719       lp->presolve_undo->fixed_rhs[j] *= scalechange[i];
720 
721     if(lp->orig_upbo[i] < lp->infinite)     /* This is the range */
722       lp->orig_upbo[i] *= scalechange[i];
723 
724     if((lp->orig_lowbo[i] != 0) && (fabs(lp->orig_lowbo[i]) < lp->infinite))
725       lp->orig_lowbo[i] *= scalechange[i];
726   }
727 
728   set_action(&lp->spx_action, ACTION_REBASE | ACTION_REINVERT | ACTION_RECOMPUTE);
729 
730   return( TRUE );
731 }
732 
scale(lprec * lp,REAL * scaledelta)733 STATIC REAL scale(lprec *lp, REAL *scaledelta)
734 {
735   int     i, j, nz, row_count, nzOF = 0;
736   REAL    *row_max, *row_min, *scalechange = NULL, absval;
737   REAL    col_max, col_min;
738   MYBOOL  rowscaled, colscaled;
739   MATrec  *mat = lp->matA;
740   REAL    *value;
741   int     *rownr;
742 
743   if(is_scaletype(lp, SCALE_NONE))
744     return(0.0);
745 
746   if(!lp->scaling_used) {
747     allocREAL(lp, &lp->scalars, lp->sum_alloc + 1, FALSE);
748     for(i = 0; i <= lp->sum; i++) {
749       lp->scalars[i] = 1;
750     }
751     lp->scaling_used = TRUE;
752   }
753 #ifdef Paranoia
754   else
755     for(i = 0; i <= lp->sum; i++) {
756       if(lp->scalars[i] == 0)
757         report(lp, SEVERE, "scale: Zero-valued scalar found at index %d\n", i);
758     }
759 #endif
760   if(scaledelta == NULL)
761     allocREAL(lp, &scalechange, lp->sum + 1, FALSE);
762   else
763     scalechange = scaledelta;
764 
765  /* Must initialize due to computation of scaling statistic - KE */
766   for(i = 0; i <= lp->sum; i++)
767     scalechange[i] = 1;
768 
769   row_count = lp->rows;
770   allocREAL(lp, &row_max, row_count + 1, TRUE);
771   allocREAL(lp, &row_min, row_count + 1, FALSE);
772 
773   /* Initialise min and max values of rows */
774   for(i = 0; i <= row_count; i++) {
775     if(is_scaletype(lp, SCALE_MEAN))
776       row_min[i] = 0;             /* Carries the count of elements */
777     else
778       row_min[i] = lp->infinite;  /* Carries the minimum element */
779   }
780 
781   /* Calculate row scaling data */
782   for(j = 1; j <= lp->columns; j++) {
783 
784     absval = lp->orig_obj[j];
785     if(absval != 0) {
786       absval = scaled_mat(lp, absval, 0, j);
787       accumulate_for_scale(lp, &row_min[0], &row_max[0], absval);
788       nzOF++;
789     }
790 
791     i = mat->col_end[j - 1];
792     value = &(COL_MAT_VALUE(i));
793     rownr = &(COL_MAT_ROWNR(i));
794     nz = mat->col_end[j];
795     for(; i < nz;
796         i++, value += matValueStep, rownr += matRowColStep) {
797       absval = scaled_mat(lp, *value, *rownr, j);
798       accumulate_for_scale(lp, &row_min[*rownr], &row_max[*rownr], absval);
799     }
800   }
801 
802   /* Calculate scale factors for rows */
803   i = 0;
804   for(; i <= lp->rows; i++) {
805     if(i == 0)
806       nz = nzOF;
807     else
808       nz = mat_rowlength(lp->matA, i);
809     absval = minmax_to_scale(lp, row_min[i], row_max[i], nzOF);
810     if(absval == 0)
811       absval = 1;
812     scalechange[i] = absval;
813   }
814 
815   FREE(row_max);
816   FREE(row_min);
817 
818   /* Row-scale the matrix (including the objective function and Lagrangean constraints) */
819   rowscaled = scale_updaterows(lp, scalechange, TRUE);
820 
821   /* Calculate column scales */
822   i = 1;
823   for(j = 1; j <= lp->columns; j++) {
824     if(is_int(lp,j) && !is_integerscaling(lp)) { /* do not scale integer columns */
825       scalechange[lp->rows + j] = 1;
826     }
827     else {
828       col_max = 0;
829       if(is_scaletype(lp, SCALE_MEAN))
830         col_min = 0;
831       else
832         col_min = lp->infinite;
833 
834       absval = lp->orig_obj[j];
835       if(absval != 0) {
836         absval = scaled_mat(lp, absval, 0, j);
837         accumulate_for_scale(lp, &col_min, &col_max, absval);
838       }
839 
840       i = mat->col_end[j - 1];
841       value = &(COL_MAT_VALUE(i));
842       rownr = &(COL_MAT_ROWNR(i));
843       nz = mat->col_end[j];
844       for(; i < nz;
845           i++, value += matValueStep, rownr += matRowColStep) {
846         absval = scaled_mat(lp, *value, *rownr, j);
847         accumulate_for_scale(lp, &col_min, &col_max, absval);
848       }
849       nz = mat_collength(lp->matA, j);
850       if(fabs(lp->orig_obj[j]) > 0)
851         nz++;
852       scalechange[lp->rows + j] = minmax_to_scale(lp, col_min, col_max, nz);
853     }
854   }
855 
856   /* ... and then column-scale the already row-scaled matrix */
857   colscaled = scale_updatecolumns(lp, &scalechange[lp->rows], TRUE);
858 
859   /* Create a geometric mean-type measure of the extent of scaling performed; */
860   /* ideally, upon successive calls to scale() the value should converge to 0 */
861   if(rowscaled || colscaled) {
862     col_max = 0;
863     for(j = 1; j <= lp->columns; j++)
864       col_max += log(scalechange[lp->rows + j]);
865     col_max = exp(col_max/lp->columns);
866 
867     i = 0;
868     col_min = 0;
869     for(; i <= lp->rows; i++)
870       col_min += log(scalechange[i]);
871     col_min = exp(col_min/row_count);
872   }
873   else {
874     col_max = 1;
875     col_min = 1;
876   }
877 
878   if(scaledelta == NULL)
879     FREE(scalechange);
880 
881   return(1 - sqrt(col_max*col_min));
882 }
883 
finalize_scaling(lprec * lp,REAL * scaledelta)884 STATIC MYBOOL finalize_scaling(lprec *lp, REAL *scaledelta)
885 {
886   int i;
887 
888   /* Check if we should equilibrate */
889   if(is_scalemode(lp, SCALE_EQUILIBRATE) && !is_scaletype(lp, SCALE_CURTISREID)) {
890     int oldmode;
891 
892     oldmode = lp->scalemode;
893     lp->scalemode = SCALE_LINEAR + SCALE_EXTREME;
894     scale(lp, scaledelta);
895     lp->scalemode = oldmode;
896   }
897 
898   /* Check if we should prevent rounding errors */
899   if(is_scalemode(lp, SCALE_POWER2)) {
900     REAL *scalars;
901     if(scaledelta == NULL)
902       scalars = lp->scalars;
903     else
904       scalars = scaledelta;
905 
906     for(i = 0; i <= lp->sum; i++)
907       scalars[i] = roundPower2(scalars[i]);
908   }
909 
910   /* Then transfer the scalars to the model's data */
911   return( scale_rows(lp, scaledelta) && scale_columns(lp, scaledelta) );
912 
913 }
914 
auto_scale(lprec * lp)915 STATIC REAL auto_scale(lprec *lp)
916 {
917   int    n = 1;
918   REAL   scalingmetric = 0, *scalenew = NULL;
919 
920   if(lp->scaling_used &&
921      ((((lp->scalemode & SCALE_DYNUPDATE) == 0)) || (lp->bb_level > 0)))
922     return( scalingmetric);
923 
924   if(lp->scalemode != SCALE_NONE) {
925 
926     /* Allocate array for incremental scaling if appropriate */
927     if((lp->solvecount > 1) && (lp->bb_level < 1) &&
928        ((lp->scalemode & SCALE_DYNUPDATE) != 0))
929       allocREAL(lp, &scalenew, lp->sum + 1, FALSE);
930 
931     if(is_scaletype(lp, SCALE_CURTISREID)) {
932       scalingmetric = scaleCR(lp, scalenew);
933     }
934     else {
935       REAL scalinglimit, scalingdelta;
936       int  count;
937 
938       /* Integer value of scalelimit holds the maximum number of iterations; default to 1 */
939       count = (int) floor(lp->scalelimit);
940       scalinglimit = lp->scalelimit;
941       if((count == 0) || (scalinglimit == 0)) {
942         if(scalinglimit > 0)
943           count = DEF_SCALINGLIMIT;  /* A non-zero convergence has been given, default to max 5 iterations */
944         else
945           count = 1;
946       }
947       else
948         scalinglimit -= count;
949 
950       /* Scale to desired relative convergence or iteration limit */
951       n = 0;
952       scalingdelta = 1.0;
953       scalingmetric = 1.0;
954       while((n < count) && (fabs(scalingdelta) > scalinglimit)) {
955         n++;
956         scalingdelta = scale(lp, scalenew);
957         scalingmetric = scalingmetric*(1+scalingdelta);
958       }
959       scalingmetric -= 1;
960     }
961   }
962 
963   /* Update the inf norm of the elements of the matrix (excluding the OF) */
964   mat_computemax(lp->matA);
965 
966   /* Check if we really have to do scaling */
967   if(lp->scaling_used && (fabs(scalingmetric) >= lp->epsprimal))
968     /* Ok, do it */
969     finalize_scaling(lp, scalenew);
970 
971   else {
972 
973     /* Otherwise reset scaling variables */
974     if(lp->scalars != NULL) {
975       FREE(lp->scalars);
976     }
977     lp->scaling_used = FALSE;
978     lp->columns_scaled = FALSE;
979   }
980   if(scalenew != NULL)
981     FREE(scalenew);
982 
983   return(scalingmetric);
984 }
985 
unscale_columns(lprec * lp)986 STATIC void unscale_columns(lprec *lp)
987 {
988   int     i, j, nz;
989   MATrec  *mat = lp->matA;
990   REAL    *value;
991   int     *rownr, *colnr;
992 
993   if(!lp->columns_scaled)
994     return;
995 
996   /* Unscale OF */
997   for(j = 1; j <= lp->columns; j++) {
998     lp->orig_obj[j] = unscaled_mat(lp, lp->orig_obj[j], 0, j);
999   }
1000 
1001   /* Unscale mat */
1002   mat_validate(mat);
1003   nz = get_nonzeros(lp);
1004   value = &(COL_MAT_VALUE(0));
1005   rownr = &(COL_MAT_ROWNR(0));
1006   colnr = &(COL_MAT_COLNR(0));
1007   for(j = 0; j < nz;
1008       j++, value += matValueStep, rownr += matRowColStep, colnr += matRowColStep) {
1009     *value = unscaled_mat(lp, *value, *rownr, *colnr);
1010   }
1011 
1012   /* Unscale bounds as well */
1013   for(i = lp->rows + 1, j = 1; i <= lp->sum; i++, j++) {
1014     lp->orig_lowbo[i] = unscaled_value(lp, lp->orig_lowbo[i], i);
1015     lp->orig_upbo[i]  = unscaled_value(lp, lp->orig_upbo[i], i);
1016     lp->sc_lobound[j]  = unscaled_value(lp, lp->sc_lobound[j], i);
1017   }
1018 
1019   for(i = lp->rows + 1; i<= lp->sum; i++)
1020     lp->scalars[i] = 1;
1021 
1022   lp->columns_scaled = FALSE;
1023   set_action(&lp->spx_action, ACTION_REBASE | ACTION_REINVERT | ACTION_RECOMPUTE);
1024 }
1025 
undoscale(lprec * lp)1026 void undoscale(lprec *lp)
1027 {
1028   int     i, j, nz;
1029   MATrec  *mat = lp->matA;
1030   REAL    *value;
1031   int     *rownr, *colnr;
1032 
1033   if(lp->scaling_used) {
1034 
1035     /* Unscale the OF */
1036     for(j = 1; j <= lp->columns; j++) {
1037       lp->orig_obj[j] = unscaled_mat(lp, lp->orig_obj[j], 0, j);
1038     }
1039 
1040     /* Unscale the matrix */
1041     mat_validate(mat);
1042     nz = get_nonzeros(lp);
1043     value = &(COL_MAT_VALUE(0));
1044     rownr = &(COL_MAT_ROWNR(0));
1045     colnr = &(COL_MAT_COLNR(0));
1046     for(j = 0; j < nz;
1047         j++, value += matValueStep, rownr += matRowColStep, colnr += matRowColStep) {
1048       *value = unscaled_mat(lp, *value, *rownr, *colnr);
1049     }
1050 
1051     /* Unscale variable bounds */
1052     for(i = lp->rows + 1, j = 1; i <= lp->sum; i++, j++) {
1053       lp->orig_lowbo[i] = unscaled_value(lp, lp->orig_lowbo[i], i);
1054       lp->orig_upbo[i]  = unscaled_value(lp, lp->orig_upbo[i], i);
1055       lp->sc_lobound[j]  = unscaled_value(lp, lp->sc_lobound[j], i);
1056     }
1057 
1058     /* Unscale the rhs, upper and lower bounds... */
1059     for(i = 0; i <= lp->rows; i++) {
1060       lp->orig_rhs[i] = unscaled_value(lp, lp->orig_rhs[i], i);
1061       j = lp->presolve_undo->var_to_orig[i];
1062       if(j != 0)
1063         lp->presolve_undo->fixed_rhs[j] = unscaled_value(lp, lp->presolve_undo->fixed_rhs[j], i);
1064       lp->orig_lowbo[i] = unscaled_value(lp, lp->orig_lowbo[i], i);
1065       lp->orig_upbo[i] = unscaled_value(lp, lp->orig_upbo[i], i);
1066     }
1067 
1068     FREE(lp->scalars);
1069     lp->scaling_used = FALSE;
1070     lp->columns_scaled = FALSE;
1071 
1072     set_action(&lp->spx_action, ACTION_REBASE | ACTION_REINVERT | ACTION_RECOMPUTE);
1073   }
1074 }
1075 
1076