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 else { /* pass is odd, compute ColScale */
356 for(col=1; col<=colMax; col++)
357 FColScale[col]*=(1 + ekekm1 / qkm1);
358 for(col=1; col<=colMax; col++)
359 FColScale[col]+=(residual_even[col] / ((REAL) ColCount[col] * qkm1) -
360 ColScalem2[col] * ekekm1 / qkm1);
361 }
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], nz); /* nz instead of nzOF KJEI 20/05/2010 */
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