1
2 /*
3 ----------------------------------------------------------------------------------
4 Crash management routines in lp_solve v5.0+
5 ----------------------------------------------------------------------------------
6 Author: Kjell Eikland
7 Contact: kjell.eikland@broadpark.no
8 License terms: LGPL.
9
10 Requires: lp_lib.h, lp_utils.h, lp_matrix.h
11
12 Release notes:
13 v1.0.0 1 April 2004 First version.
14 v1.1.0 20 July 2004 Reworked with flexible matrix storage model.
15
16 ----------------------------------------------------------------------------------
17 */
18
19 #include <string.h>
20
21 #include "commonlib.h"
22 #include "lp_lib.h"
23 #include "lp_scale.h"
24 #include "lp_utils.h"
25 #include "lp_report.h"
26 #include "lp_matrix.h"
27 #include "lp_crash.h"
28
29 #ifdef FORTIFY
30 # include "lp_fortify.h"
31 #endif
32
33
crash_basis(lprec * lp)34 MYBOOL crash_basis(lprec *lp)
35 {
36 int i;
37 MATrec *mat = lp->matA;
38 MYBOOL ok = TRUE;
39
40 /* Initialize basis indicators */
41 if(lp->basis_valid)
42 lp->var_basic[0] = FALSE;
43 else
44 default_basis(lp);
45
46 /* Set initial partial pricing blocks */
47 if(lp->rowblocks != NULL)
48 lp->rowblocks->blocknow = 1;
49 if(lp->colblocks != NULL)
50 lp->colblocks->blocknow = ((lp->crashmode == CRASH_NONE) || (lp->colblocks->blockcount == 1) ? 1 : 2);
51
52 /* Construct a basis that is in some measure the "most feasible" */
53 if((lp->crashmode == CRASH_MOSTFEASIBLE) && mat_validate(mat)) {
54 /* The logic here follows Maros */
55 LLrec *rowLL = NULL, *colLL = NULL;
56 int ii, rx, cx, ix, nz;
57 REAL wx, tx, *rowMAX = NULL, *colMAX = NULL;
58 int *rowNZ = NULL, *colNZ = NULL, *rowWT = NULL, *colWT = NULL;
59 REAL *value;
60 int *rownr, *colnr;
61
62 report(lp, NORMAL, "crash_basis: 'Most feasible' basis crashing selected\n");
63
64 /* Tally row and column non-zero counts */
65 ok = allocINT(lp, &rowNZ, lp->rows+1, TRUE) &&
66 allocINT(lp, &colNZ, lp->columns+1, TRUE) &&
67 allocREAL(lp, &rowMAX, lp->rows+1, FALSE) &&
68 allocREAL(lp, &colMAX, lp->columns+1, FALSE);
69 if(!ok)
70 goto Finish;
71
72 nz = mat_nonzeros(mat);
73 rownr = &COL_MAT_ROWNR(0);
74 colnr = &COL_MAT_COLNR(0);
75 value = &COL_MAT_VALUE(0);
76 for(i = 0; i < nz;
77 i++, rownr += matRowColStep, colnr += matRowColStep, value += matValueStep) {
78 rx = *rownr;
79 cx = *colnr;
80 wx = fabs(*value);
81 rowNZ[rx]++;
82 colNZ[cx]++;
83 if(i == 0) {
84 rowMAX[rx] = wx;
85 colMAX[cx] = wx;
86 colMAX[0] = wx;
87 }
88 else {
89 SETMAX(rowMAX[rx], wx);
90 SETMAX(colMAX[cx], wx);
91 SETMAX(colMAX[0], wx);
92 }
93 }
94 /* Reduce counts for small magnitude to preserve stability */
95 rownr = &COL_MAT_ROWNR(0);
96 colnr = &COL_MAT_COLNR(0);
97 value = &COL_MAT_VALUE(0);
98 for(i = 0; i < nz;
99 i++, rownr += matRowColStep, colnr += matRowColStep, value += matValueStep) {
100 rx = *rownr;
101 cx = *colnr;
102 wx = fabs(*value);
103 #ifdef CRASH_SIMPLESCALE
104 if(wx < CRASH_THRESHOLD * colMAX[0]) {
105 rowNZ[rx]--;
106 colNZ[cx]--;
107 }
108 #else
109 if(wx < CRASH_THRESHOLD * rowMAX[rx])
110 rowNZ[rx]--;
111 if(wx < CRASH_THRESHOLD * colMAX[cx])
112 colNZ[cx]--;
113 #endif
114 }
115
116 /* Set up priority tables */
117 ok = allocINT(lp, &rowWT, lp->rows+1, TRUE);
118 createLink(lp->rows, &rowLL, NULL);
119 ok &= (rowLL != NULL);
120 if(!ok)
121 goto Finish;
122 for(i = 1; i <= lp->rows; i++) {
123 if(get_constr_type(lp, i)==EQ)
124 ii = 3;
125 else if(lp->upbo[i] < lp->infinite)
126 ii = 2;
127 else if(fabs(lp->rhs[i]) < lp->infinite)
128 ii = 1;
129 else
130 ii = 0;
131 rowWT[i] = ii;
132 if(ii > 0)
133 appendLink(rowLL, i);
134 }
135 ok = allocINT(lp, &colWT, lp->columns+1, TRUE);
136 createLink(lp->columns, &colLL, NULL);
137 ok &= (colLL != NULL);
138 if(!ok)
139 goto Finish;
140 for(i = 1; i <= lp->columns; i++) {
141 ix = lp->rows+i;
142 if(is_unbounded(lp, i))
143 ii = 3;
144 else if(lp->upbo[ix] >= lp->infinite)
145 ii = 2;
146 else if(fabs(lp->upbo[ix]-lp->lowbo[ix]) > lp->epsmachine)
147 ii = 1;
148 else
149 ii = 0;
150 colWT[i] = ii;
151 if(ii > 0)
152 appendLink(colLL, i);
153 }
154
155 /* Loop over all basis variables */
156 for(i = 1; i <= lp->rows; i++) {
157
158 /* Select row */
159 rx = 0;
160 wx = -lp->infinite;
161 for(ii = firstActiveLink(rowLL); ii > 0; ii = nextActiveLink(rowLL, ii)) {
162 tx = rowWT[ii] - CRASH_SPACER*rowNZ[ii];
163 if(tx > wx) {
164 rx = ii;
165 wx = tx;
166 }
167 }
168 if(rx == 0)
169 break;
170 removeLink(rowLL, rx);
171
172 /* Select column */
173 cx = 0;
174 wx = -lp->infinite;
175 for(ii = mat->row_end[rx-1]; ii < mat->row_end[rx]; ii++) {
176
177 /* Update NZ column counts for row selected above */
178 tx = fabs(ROW_MAT_VALUE(ii));
179 ix = ROW_MAT_COLNR(ii);
180 #ifdef CRASH_SIMPLESCALE
181 if(tx >= CRASH_THRESHOLD * colMAX[0])
182 #else
183 if(tx >= CRASH_THRESHOLD * colMAX[ix])
184 #endif
185 colNZ[ix]--;
186 if(!isActiveLink(colLL, ix) || (tx < CRASH_THRESHOLD * rowMAX[rx]))
187 continue;
188
189 /* Now do the test for best pivot */
190 tx = my_sign(lp->orig_obj[ix]) - my_sign(ROW_MAT_VALUE(ii));
191 tx = colWT[ix] + CRASH_WEIGHT*tx - CRASH_SPACER*colNZ[ix];
192 if(tx > wx) {
193 cx = ix;
194 wx = tx;
195 }
196 }
197 if(cx == 0)
198 break;
199 removeLink(colLL, cx);
200
201 /* Update row NZ counts */
202 ii = mat->col_end[cx-1];
203 rownr = &COL_MAT_ROWNR(ii);
204 value = &COL_MAT_VALUE(ii);
205 for(; ii < mat->col_end[cx];
206 ii++, rownr += matRowColStep, value += matValueStep) {
207 wx = fabs(*value);
208 ix = *rownr;
209 #ifdef CRASH_SIMPLESCALE
210 if(wx >= CRASH_THRESHOLD * colMAX[0])
211 #else
212 if(wx >= CRASH_THRESHOLD * rowMAX[ix])
213 #endif
214 rowNZ[ix]--;
215 }
216
217 /* Set new basis variable */
218 set_basisvar(lp, rx, lp->rows+cx);
219 }
220
221 /* Clean up */
222 Finish:
223 FREE(rowNZ);
224 FREE(colNZ);
225 FREE(rowMAX);
226 FREE(colMAX);
227 FREE(rowWT);
228 FREE(colWT);
229 freeLink(&rowLL);
230 freeLink(&colLL);
231 }
232
233 /* Construct a basis that is in some measure the "least degenerate" */
234 else if((lp->crashmode == CRASH_LEASTDEGENERATE) && mat_validate(mat)) {
235 /* The logic here follows Maros */
236 LLrec *rowLL = NULL, *colLL = NULL;
237 int ii, rx, cx, ix, nz, *merit = NULL;
238 REAL *value, wx, hold, *rhs = NULL, *eta = NULL;
239 int *rownr, *colnr;
240
241 report(lp, NORMAL, "crash_basis: 'Least degenerate' basis crashing selected\n");
242
243 /* Create temporary arrays */
244 ok = allocINT(lp, &merit, lp->columns + 1, FALSE) &&
245 allocREAL(lp, &eta, lp->rows + 1, FALSE) &&
246 allocREAL(lp, &rhs, lp->rows + 1, FALSE);
247 createLink(lp->columns, &colLL, NULL);
248 createLink(lp->rows, &rowLL, NULL);
249 ok &= (colLL != NULL) && (rowLL != NULL);
250 if(!ok)
251 goto FinishLD;
252 MEMCOPY(rhs, lp->orig_rhs, lp->rows + 1);
253 for(i = 1; i <= lp->columns; i++)
254 appendLink(colLL, i);
255 for(i = 1; i <= lp->rows; i++)
256 appendLink(rowLL, i);
257
258 /* Loop until we have found enough new bases */
259 while(colLL->count > 0) {
260
261 /* Tally non-zeros matching in RHS and each active column */
262 nz = mat_nonzeros(mat);
263 rownr = &COL_MAT_ROWNR(0);
264 colnr = &COL_MAT_COLNR(0);
265 ii = 0;
266 MEMCLEAR(merit, lp->columns + 1);
267 for(i = 0; i < nz;
268 i++, rownr += matRowColStep, colnr += matRowColStep) {
269 rx = *rownr;
270 cx = *colnr;
271 if(isActiveLink(colLL, cx) && (rhs[rx] != 0)) {
272 merit[cx]++;
273 ii++;
274 }
275 }
276 if(ii == 0)
277 break;
278
279 /* Find maximal match; break ties with column length */
280 i = firstActiveLink(colLL);
281 cx = i;
282 for(i = nextActiveLink(colLL, i); i != 0; i = nextActiveLink(colLL, i)) {
283 if(merit[i] >= merit[cx]) {
284 if((merit[i] > merit[cx]) || (mat_collength(mat, i) > mat_collength(mat, cx)))
285 cx = i;
286 }
287 }
288
289 /* Determine the best pivot row */
290 i = mat->col_end[cx-1];
291 nz = mat->col_end[cx];
292 rownr = &COL_MAT_ROWNR(i);
293 value = &COL_MAT_VALUE(i);
294 rx = 0;
295 wx = 0;
296 MEMCLEAR(eta, lp->rows + 1);
297 for(; i < nz;
298 i++, rownr += matRowColStep, value += matValueStep) {
299 ix = *rownr;
300 hold = *value;
301 eta[ix] = rhs[ix] / hold;
302 hold = fabs(hold);
303 if(isActiveLink(rowLL, ix) && (hold > wx)) {
304 wx = hold;
305 rx = ix;
306 }
307 }
308
309 /* Set new basis variable */
310 if(rx > 0) {
311
312 /* We have to update the rhs vector for the implied transformation
313 in order to be able to find the new RHS non-zero pattern */
314 for(i = 1; i <= lp->rows; i++)
315 rhs[i] -= wx * eta[i];
316 rhs[rx] = wx;
317
318 /* Do the exchange */
319 set_basisvar(lp, rx, lp->rows+cx);
320 removeLink(rowLL, rx);
321 }
322 removeLink(colLL, cx);
323
324 }
325
326 /* Clean up */
327 FinishLD:
328 FREE(merit);
329 FREE(rhs);
330 freeLink(&rowLL);
331 freeLink(&colLL);
332
333 }
334 return( ok );
335 }
336
337 #if 0
338 MYBOOL __WINAPI guess_basis(lprec *lp, REAL *guessvector, int *basisvector)
339 {
340 MYBOOL status = FALSE;
341 REAL *values = NULL, *violation = NULL,
342 *value, error, upB, loB, sortorder = 1.0;
343 int i, n, *rownr, *colnr;
344 MATrec *mat = lp->matA;
345
346 if(!mat_validate(lp->matA))
347 return( status );
348
349 /* Create helper arrays */
350 if(!allocREAL(lp, &values, lp->sum+1, TRUE) ||
351 !allocREAL(lp, &violation, lp->sum+1, TRUE))
352 goto Finish;
353
354 /* Compute values of slack variables for given guess vector */
355 i = 0;
356 n = get_nonzeros(lp);
357 rownr = &COL_MAT_ROWNR(i);
358 colnr = &COL_MAT_COLNR(i);
359 value = &COL_MAT_VALUE(i);
360 for(; i < n; i++, rownr += matRowColStep, colnr += matRowColStep, value += matValueStep)
361 values[*rownr] += unscaled_mat(lp, my_chsign(is_chsign(lp, *rownr), *value), *rownr, *colnr) *
362 guessvector[*colnr];
363 MEMMOVE(values+lp->rows+1, guessvector+1, lp->columns);
364
365 /* Initialize constraint bound violation measures */
366 for(i = 1; i <= lp->rows; i++) {
367 upB = get_rh_upper(lp, i);
368 loB = get_rh_lower(lp, i);
369 error = values[i] - upB;
370 if(error > lp->epsprimal)
371 violation[i] = sortorder*error;
372 else {
373 error = loB - values[i];
374 if(error > lp->epsprimal)
375 violation[i] = sortorder*error;
376 else if(is_infinite(lp, loB) && is_infinite(lp, upB))
377 ;
378 else if(is_infinite(lp, upB))
379 violation[i] = sortorder*(loB - values[i]);
380 else if(is_infinite(lp, loB))
381 violation[i] = sortorder*(values[i] - upB);
382 else
383 violation[i] = - sortorder*MAX(upB - values[i], values[i] - loB);
384 }
385 basisvector[i] = i;
386 }
387
388 /* Initialize user variable bound violation measures */
389 for(i = 1; i <= lp->columns; i++) {
390 n = lp->rows+i;
391 upB = get_upbo(lp, i);
392 loB = get_lowbo(lp, i);
393 error = guessvector[i] - upB;
394 if(error > lp->epsprimal)
395 violation[n] = sortorder*error;
396 else {
397 error = loB - values[n];
398 if(error > lp->epsprimal)
399 violation[n] = sortorder*error;
400 else if(is_infinite(lp, loB) && is_infinite(lp, upB))
401 ;
402 else if(is_infinite(lp, upB))
403 violation[n] = sortorder*(loB - values[n]);
404 else if(is_infinite(lp, loB))
405 violation[n] = sortorder*(values[n] - upB);
406 else
407 violation[n] = - sortorder*MAX(upB - values[n], values[n] - loB);
408 }
409 basisvector[n] = n;
410 }
411
412 /* Sort decending by violation; this means that variables with
413 the largest violations will be designated as basic */
414 sortByREAL(basisvector, violation, lp->sum, 1, FALSE);
415
416 /* Adjust the non-basic indeces for the (proximal) bound state */
417 error = lp->epsprimal;
418 for(i = lp->rows+1, rownr = basisvector+i; i <= lp->sum; i++, rownr++) {
419 if(*rownr <= lp->rows) {
420 if(values[*rownr] <= get_rh_lower(lp, *rownr)+error)
421 *rownr = -(*rownr);
422 }
423 else
424 if(values[i] <= get_lowbo(lp, (*rownr)-lp->rows)+error)
425 *rownr = -(*rownr);
426 }
427
428 /* Clean up and return status */
429 status = (MYBOOL) (violation[1] == 0);
430 Finish:
431 FREE(values);
432 FREE(violation);
433
434
435 return( status );
436 }
437 #endif
438
439 #if 0
440 MYBOOL __WINAPI guess_basis(lprec *lp, REAL *guessvector, int *basisvector)
441 {
442 MYBOOL *isnz, status = FALSE;
443 REAL *values = NULL, *violation = NULL,
444 eps = lp->epsprimal,
445 *value, error, upB, loB, sortorder = 1.0;
446 int i, j, n, *rownr, *colnr, *slkpos,
447 nrows = lp->rows, ncols = lp->columns;
448 MATrec *mat = lp->matA;
449
450 if(!mat_validate(mat))
451 return( status );
452
453 /* Create helper arrays */
454 if(!allocREAL(lp, &values, lp->sum+1, TRUE) ||
455 !allocREAL(lp, &violation, lp->sum+1, TRUE))
456 goto Finish;
457
458 /* Compute values of slack variables for given guess vector */
459 i = 0;
460 n = get_nonzeros(lp);
461 rownr = &COL_MAT_ROWNR(i);
462 colnr = &COL_MAT_COLNR(i);
463 value = &COL_MAT_VALUE(i);
464 for(; i < n; i++, rownr += matRowColStep, colnr += matRowColStep, value += matValueStep)
465 values[*rownr] += unscaled_mat(lp, my_chsign(is_chsign(lp, *rownr), *value), *rownr, *colnr) *
466 guessvector[*colnr];
467 MEMMOVE(values+nrows+1, guessvector+1, ncols);
468
469 /* Initialize constraint bound violation measures (expressed as positive values) */
470 for(i = 1; i <= nrows; i++) {
471 upB = get_rh_upper(lp, i);
472 loB = get_rh_lower(lp, i);
473 error = values[i] - upB;
474 if(error > eps)
475 violation[i] = sortorder*error;
476 else {
477 error = loB - values[i];
478 if(error > eps)
479 violation[i] = sortorder*error;
480 else if(my_infinite(lp, loB) && my_infinite(lp, upB))
481 ;
482 else if(my_infinite(lp, upB))
483 violation[i] = sortorder*(loB - values[i]);
484 else if(my_infinite(lp, loB))
485 violation[i] = sortorder*(values[i] - upB);
486 else
487 violation[i] = -sortorder*MAX(upB - values[i], values[i] - loB);
488 }
489 basisvector[i] = i;
490 }
491
492 /* Initialize user variable bound violation measures (expressed as positive values) */
493 for(i = 1; i <= ncols; i++) {
494 n = nrows+i;
495 upB = get_upbo(lp, i);
496 loB = get_lowbo(lp, i);
497 error = guessvector[i] - upB;
498 if(error > eps)
499 violation[n] = sortorder*error;
500 else {
501 error = loB - values[n];
502 if(error > eps)
503 violation[n] = sortorder*error;
504 else if(my_infinite(lp, loB) && my_infinite(lp, upB))
505 ;
506 else if(my_infinite(lp, upB))
507 violation[n] = sortorder*(loB - values[n]);
508 else if(my_infinite(lp, loB))
509 violation[n] = sortorder*(values[n] - upB);
510 else
511 violation[n] = -sortorder*MAX(upB - values[n], values[n] - loB);
512 }
513 basisvector[n] = n;
514 }
515
516 /* Sort decending by violation; this means that variables with
517 the largest violations will be designated as basic */
518 sortByREAL(basisvector, violation, lp->sum, 1, FALSE);
519 error = violation[1];
520
521 /* Adjust the non-basic indeces for the (proximal) bound state */
522 for(i = nrows+1, rownr = basisvector+i; i <= lp->sum; i++, rownr++) {
523 if(*rownr <= nrows) {
524 if(values[*rownr] <= get_rh_lower(lp, *rownr)+eps)
525 *rownr = -(*rownr);
526 }
527 else
528 if(values[i] <= get_lowbo(lp, (*rownr)-nrows)+eps)
529 *rownr = -(*rownr);
530 }
531
532 #if 1
533 /* Let us check for obvious row singularities and try to fix these;
534 First assemble necessary basis statistics... */
535 isnz = (MYBOOL *) values;
536 MEMCLEAR(isnz, nrows+1);
537 slkpos = (int *) violation;
538 MEMCLEAR(slkpos, nrows+1);
539 for(i = 1; i <= nrows; i++) {
540 j = abs(basisvector[i]);
541 if(j <= nrows) {
542 isnz[j] = TRUE;
543 slkpos[j] = i;
544 }
545 else {
546 j-= nrows;
547 j = mat->col_end[j-1];
548 isnz[COL_MAT_ROWNR(j)] = TRUE;
549 /*isnz[COL_MAT_ROWNR(j+1)] = TRUE;*/
550 }
551 }
552 for(; i <= lp->sum; i++) {
553 j = abs(basisvector[i]);
554 if(j <= nrows)
555 slkpos[j] = i;
556 }
557
558 /* ...then set the corresponding slacks basic for row rank deficient positions */
559 for(j = 1; j <= nrows; j++) {
560 #ifdef Paranoia
561 if(slkpos[j] == 0)
562 report(lp, SEVERE, "guess_basis: Internal error");
563 #endif
564 if(!isnz[j]) {
565 isnz[j] = TRUE;
566 i = slkpos[j];
567 swapINT(&basisvector[i], &basisvector[j]);
568 basisvector[j] = abs(basisvector[j]);
569 }
570 }
571 #endif
572
573 /* Clean up and return status */
574 status = (MYBOOL) (error <= eps);
575 Finish:
576 FREE(values);
577 FREE(violation);
578
579 return( status );
580 }
581 #endif
582
583 #if 0
584 MYBOOL __WINAPI guess_basis(lprec *lp, REAL *guessvector, int *basisvector)
585 {
586 MYBOOL *isnz, status = FALSE;
587 REAL *values = NULL, *violation = NULL,
588 eps = lp->epsprimal,
589 *value, error, upB, loB, sortorder = 1.0;
590 int i, j, jj, n, *rownr, *colnr, *slkpos,
591 nrows = lp->rows, ncols = lp->columns;
592 MATrec *mat = lp->matA;
593
594 if(!mat_validate(mat))
595 return( status );
596
597 /* Create helper arrays */
598 if(!allocREAL(lp, &values, lp->sum+1, TRUE) ||
599 !allocREAL(lp, &violation, lp->sum+1, TRUE))
600 goto Finish;
601
602 /* Compute values of slack variables for given guess vector */
603 i = 0;
604 n = get_nonzeros(lp);
605 rownr = &COL_MAT_ROWNR(i);
606 colnr = &COL_MAT_COLNR(i);
607 value = &COL_MAT_VALUE(i);
608 for(; i < n; i++, rownr += matRowColStep, colnr += matRowColStep, value += matValueStep)
609 values[*rownr] += unscaled_mat(lp, my_chsign(is_chsign(lp, *rownr), *value), *rownr, *colnr) *
610 guessvector[*colnr];
611 MEMMOVE(values+nrows+1, guessvector+1, ncols);
612
613 /* Initialize constraint bound violation measures (expressed as positive values) */
614 for(i = 1; i <= nrows; i++) {
615 upB = get_rh_upper(lp, i);
616 loB = get_rh_lower(lp, i);
617 error = values[i] - upB;
618 if(error > -eps)
619 violation[i] = sortorder*MAX(0,error);
620 else {
621 error = loB - values[i];
622 if(error > -eps)
623 violation[i] = sortorder*MAX(0,error);
624 else if(my_infinite(lp, loB) && my_infinite(lp, upB))
625 ;
626 else if(my_infinite(lp, upB))
627 violation[i] = sortorder*(loB - values[i]);
628 else if(my_infinite(lp, loB))
629 violation[i] = sortorder*(values[i] - upB);
630 else
631 violation[i] = -sortorder*MAX(upB - values[i], values[i] - loB);
632 }
633 basisvector[i] = i;
634 }
635
636 /* Initialize user variable bound violation measures (expressed as positive values) */
637 for(i = 1; i <= ncols; i++) {
638 n = nrows+i;
639 upB = get_upbo(lp, i);
640 loB = get_lowbo(lp, i);
641 error = guessvector[i] - upB;
642 if(error > -eps)
643 violation[n] = sortorder*MAX(0,error);
644 else {
645 error = loB - values[n];
646 if(error > -eps)
647 violation[n] = sortorder*MAX(0,error);
648 else if(my_infinite(lp, loB) && my_infinite(lp, upB))
649 ;
650 else if(my_infinite(lp, upB))
651 violation[n] = sortorder*(loB - values[n]);
652 else if(my_infinite(lp, loB))
653 violation[n] = sortorder*(values[n] - upB);
654 else
655 violation[n] = -sortorder*MAX(upB - values[n], values[n] - loB);
656 }
657 basisvector[n] = n;
658 }
659
660 /* Sort decending by violation; this means that variables with
661 the largest violations will be designated as basic */
662 sortByREAL(basisvector, violation, lp->sum, 1, FALSE);
663 error = violation[1];
664
665 /* Adjust the non-basic indeces for the (proximal) bound state */
666 for(i = nrows+1, rownr = basisvector+i; i <= lp->sum; i++, rownr++) {
667 if(*rownr <= nrows) {
668 values[*rownr] -= lp->orig_rhs[*rownr];
669 if(values[*rownr] <= eps)
670 *rownr = -(*rownr);
671 }
672 else
673 if(values[i] <= get_lowbo(lp, (*rownr)-nrows)+eps)
674 *rownr = -(*rownr);
675 }
676
677 /* Let us check for obvious row singularities and try to fix these;
678 First assemble necessary basis statistics... */
679 isnz = (MYBOOL *) values;
680 MEMCLEAR(isnz, nrows+1);
681 slkpos = (int *) violation;
682 MEMCLEAR(slkpos, nrows+1);
683 for(i = 1; i <= nrows; i++) {
684 j = abs(basisvector[i]);
685 if(j <= nrows) {
686 isnz[j] = TRUE;
687 slkpos[j] = i;
688 }
689 else {
690 j-= nrows;
691 jj = mat->col_end[j-1];
692 isnz[COL_MAT_ROWNR(jj)] = TRUE;
693 /* if(++jj < mat->col_end[j])
694 isnz[COL_MAT_ROWNR(jj)] = TRUE; */
695 }
696 }
697 for(; i <= lp->sum; i++) {
698 j = abs(basisvector[i]);
699 if(j <= nrows)
700 slkpos[j] = i;
701 }
702
703 /* ...then set the corresponding slacks basic for row rank deficient positions */
704 for(j = 1; j <= nrows; j++) {
705 #ifdef Paranoia
706 if(slkpos[j] == 0)
707 report(lp, SEVERE, "guess_basis: Internal error");
708 #endif
709 if(!isnz[j]) {
710 isnz[j] = TRUE;
711 i = slkpos[j];
712 swapINT(&basisvector[i], &basisvector[j]);
713 basisvector[j] = abs(basisvector[j]);
714 }
715 }
716
717 /* Lastly normalize all basic variables to be coded as lower-bounded */
718 for(i = 1; i <= nrows; i++)
719 basisvector[i] = -abs(basisvector[i]);
720
721 /* Clean up and return status */
722 status = (MYBOOL) (error <= eps);
723 Finish:
724 FREE(values);
725 FREE(violation);
726
727 return( status );
728 }
729 #endif
730
guess_basis(lprec * lp,REAL * guessvector,int * basisvector)731 MYBOOL __WINAPI guess_basis(lprec *lp, REAL *guessvector, int *basisvector)
732 {
733 MYBOOL *isnz = NULL, status = FALSE;
734 REAL *values = NULL, *violation = NULL,
735 eps = lp->epsprimal,
736 *value, error, upB, loB, sortorder = -1.0;
737 int i, j, jj, n, *rownr, *colnr, *slkpos = NULL,
738 nrows = lp->rows, ncols = lp->columns, nsum = lp->sum;
739 int *basisnr;
740 MATrec *mat = lp->matA;
741
742 if(!mat_validate(mat))
743 return( status );
744
745 /* Create helper arrays, providing for multiple use of the violation array */
746 if(!allocREAL(lp, &values, nsum+1, TRUE) ||
747 !allocREAL(lp, &violation, nsum+1, TRUE))
748 goto Finish;
749
750 /* Compute the values of the constraints for the given guess vector */
751 i = 0;
752 n = get_nonzeros(lp);
753 rownr = &COL_MAT_ROWNR(i);
754 colnr = &COL_MAT_COLNR(i);
755 value = &COL_MAT_VALUE(i);
756 for(; i < n; i++, rownr += matRowColStep, colnr += matRowColStep, value += matValueStep)
757 values[*rownr] += unscaled_mat(lp, my_chsign(is_chsign(lp, *rownr), *value), *rownr, *colnr) *
758 guessvector[*colnr];
759 MEMMOVE(values+nrows+1, guessvector+1, ncols);
760
761 /* Initialize bound "violation" or primal non-degeneracy measures, expressed
762 as the absolute value of the differences from the closest bound. */
763 for(i = 1; i <= nsum; i++) {
764 if(i <= nrows) {
765 loB = get_rh_lower(lp, i);
766 upB = get_rh_upper(lp, i);
767 }
768 else {
769 loB = get_lowbo(lp, i-nrows);
770 upB = get_upbo(lp, i-nrows);
771 }
772
773 /* Free constraints/variables */
774 if(my_infinite(lp, loB) && my_infinite(lp, upB))
775 error = 0;
776 /* Violated constraints/variable bounds */
777 else if(values[i]+eps < loB)
778 error = loB-values[i];
779 else if(values[i]-eps > upB)
780 error = values[i]-upB;
781 /* Non-violated constraints/variables bounds */
782 else if(my_infinite(lp, upB))
783 error = MAX(0, values[i]-loB);
784 else if(my_infinite(lp, loB))
785 error = MAX(0, upB-values[i]);
786 else
787 error = MIN(upB-values[i], values[i]-loB); /* MAX(upB-values[i], values[i]-loB); */
788 if(error != 0)
789 violation[i] = sortorder*error;
790 basisvector[i] = i;
791 }
792
793 /* Sort decending , meaning that variables with the largest
794 "violations" will be designated basic. Effectively, we are performing a
795 greedy type algorithm, but start at the "least interesting" end. */
796 sortByREAL(basisvector, violation, nsum, 1, FALSE);
797 error = violation[1]; /* Used for setting the return value */
798
799 /* Let us check for obvious row singularities and try to fix these.
800 Note that we reuse the memory allocated to the violation array.
801 First assemble necessary basis statistics... */
802 slkpos = (int *) violation;
803 n = nrows+1;
804 MEMCLEAR(slkpos, n);
805 isnz = (MYBOOL *) (slkpos+n+1);
806 MEMCLEAR(isnz, n);
807 for(i = 1; i <= nrows; i++) {
808 j = abs(basisvector[i]);
809 if(j <= nrows) {
810 isnz[j] = TRUE;
811 slkpos[j] = i;
812 }
813 else {
814 j-= nrows;
815 jj = mat->col_end[j-1];
816 jj = COL_MAT_ROWNR(jj);
817 isnz[jj] = TRUE;
818 }
819 }
820 for(; i <= nsum; i++) {
821 j = abs(basisvector[i]);
822 if(j <= nrows)
823 slkpos[j] = i;
824 }
825
826 /* ...then set the corresponding slacks basic for row rank deficient positions */
827 for(j = 1; j <= nrows; j++) {
828 if(slkpos[j] == 0)
829 report(lp, SEVERE, "guess_basis: Internal error");
830 if(!isnz[j]) {
831 isnz[j] = TRUE;
832 i = slkpos[j];
833 swapINT(&basisvector[i], &basisvector[j]);
834 basisvector[j] = abs(basisvector[j]);
835 }
836 }
837
838 /* Adjust the non-basic indeces for the (proximal) bound state */
839 for(i = nrows+1, basisnr = basisvector+i; i <= nsum; i++, basisnr++) {
840 n = *basisnr;
841 if(n <= nrows) {
842 values[n] -= get_rh_lower(lp, n);
843 if(values[n] <= eps)
844 *basisnr = -(*basisnr);
845 }
846 else
847 if(values[n]-eps <= get_lowbo(lp, n-nrows))
848 *basisnr = -(*basisnr);
849 }
850
851 /* Lastly normalize all basic variables to be coded as lower-bounded,
852 or effectively zero-based in the case of free variables. */
853 for(i = 1; i <= nrows; i++)
854 basisvector[i] = -abs(basisvector[i]);
855
856 /* Clean up and return status */
857 status = (MYBOOL) (error <= eps);
858 Finish:
859 FREE(values);
860 FREE(violation);
861
862 return( status );
863 }
864