1/* $Id: CoinPresolvePsdebug.cpp 1373 2011-01-03 23:57:44Z lou $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others. All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include <new>
7#include <stdio.h>
8#include <math.h>
9
10#include "CoinPresolveMatrix.hpp"
11#include "CoinHelperFunctions.hpp"
12
13/*
14 \file
15
16 This file contains a number of routines that are useful when doing serious
17 debugging but unneeded otherwise. Presumably if you're deep enough into
18 presolve to need these, you're willing to scan the file to see what can be
19 done. See also the Presolve Debug Functions module in the doxygen doc'n
20 and CoinPresolvePsdebug.hpp.
21
22 The general approach for the matrix consistency routines is that the
23 routines return void and abort when they find a problem. The routines
24 that check the basis and solution complain loudly but do not abort.
25
26 NOTE: The definitions for PRESOLVE_CONSISTENCY and PRESOLVE_DEBUG MUST BE
27 CONSISTENT across all CoinPresolve source files. AND OsiPresolve,
28 if you're debugging there. Otherwise, at best you'll get garbage
29 output. More likely, you'll get a core dump. Resist the temptation to
30 define these constants in individual files. In particular, cdone and
31 rdone will NOT be consistently maintained during postsolve.
32
33 Hack away as your needs dictate.
34*/
35
36
37/*
38 Integrity checking routines for the (loosely) packed matrices of a
39 CoinPresolveMatrix object. Some routines work on column-major and row-major
40 reps separately, others do cross-checking.
41*/
42
43namespace { // begin unnamed file-local namespace
44
45#if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
46/*
47 Check for duplicate entries in a major vector by walking the vector. For each
48 coefficient, use presolve_find_row to search the remainder of the column
49 for an entry with the same row index. We don't want to find anything.
50*/
51
52void no_majvec_dups (const char *majdones, const CoinBigIndex *majstrts,
53 const int *minndxs, const int *majlens, int nmaj)
54
55{ for (int maj = 0 ; maj < nmaj ; maj++)
56 { if ((!majdones || majdones[maj]) && majlens[maj] > 0)
57 { CoinBigIndex ks = majstrts[maj] ;
58 CoinBigIndex ke = ks+majlens[maj] ;
59 for (CoinBigIndex k = ks ; k < ke ; k++)
60 {
61/*
62 Assert we fell off the end of the column without finding the entry.
63*/
64 PRESOLVEASSERT(presolve_find_minor1(minndxs[k],k+1,
65 ke,minndxs) == ke) ; } } }
66 return ; }
67
68/*
69 As the name implies: scan for explicit zeros.
70*/
71void check_majvec_nozeros (const CoinBigIndex *majstrts, const double *majels,
72 const int *majlens, int nmaj)
73
74{ for (int maj = 0 ; maj < nmaj ; maj++)
75 { if (majlens[maj] > 0)
76 { CoinBigIndex ks = majstrts[maj] ;
77 CoinBigIndex ke = ks+majlens[maj] ;
78 for (CoinBigIndex k = ks ; k < ke ; k++)
79 { PRESOLVEASSERT(fabs(majels[k]) > ZTOLDP) ; } } }
80
81 return ; }
82
83/*
84 Integrity checks for the linked lists that indicate major vector ordering
85 in the bulk storage area (minor index and coefficient arrays).
86 */
87void links_ok (presolvehlink *majlink, int *majstrts, int *majlens, int nmaj)
88
89{ int maj ;
90
91/*
92 Confirm link integrity. Vectors of length 0 should not be part of the chain.
93*/
94 for (maj = 0 ; maj < nmaj ; maj++)
95 { int pre = majlink[maj].pre ;
96 int suc = majlink[maj].suc ;
97
98 if (majlens[maj] == 0)
99 { PRESOLVEASSERT(pre == NO_LINK && suc == NO_LINK) ; }
100 if (pre != NO_LINK)
101 { PRESOLVEASSERT(0 <= pre && pre <= nmaj) ;
102 PRESOLVEASSERT(majlink[pre].suc == maj) ; }
103 if (suc != NO_LINK)
104 { PRESOLVEASSERT(0 <= suc && suc <= nmaj) ;
105 PRESOLVEASSERT(majlink[suc].pre == maj) ; } }
106/*
107 There must be a first vector.
108*/
109 for (maj = 0 ; maj < nmaj ; maj++)
110 { if (majlink[maj].pre == NO_LINK)
111 break ; }
112 PRESOLVEASSERT(nmaj == 0 || maj < nmaj) ;
113/*
114 The order of the linked list should match the ordering indicated by the
115 major vector start & length arrays.
116*/
117 while (maj != NO_LINK)
118 { if (majlink[maj].suc != NO_LINK)
119 { PRESOLVEASSERT(majstrts[maj]+majlens[maj] <=
120 majstrts[majlink[maj].suc]) ; }
121 maj = majlink[maj].suc ; }
122
123 return ; }
124
125
126/*
127 matrix_consistent checks that an entry is in the column-major representation
128 if it is in the row-major representation. If testvals is non-zero, it also
129 checks that their values are the same.
130
131 By doing the appropriate swaps of column- and row-major data structures in
132 the parameter list, we can check that an entry is in the row-major
133 representation if it's in the column-major representation.
134
135 I can't see any nice way to rename the parameters (majmajstrt? minmajstrt?).
136
137 Original comment: ``Note that there may be entries in a row that correspond
138 to empty columns and vice-versa.'' To which a previous browser had commented
139 ``HUH???''. And I agree. -- lh, 040907 --
140*/
141
142void matrix_consistent (const CoinBigIndex *mrstrt, const int *hinrow,
143 const int *hcol, const double *rowels,
144 const CoinBigIndex *mcstrt, const int *hincol,
145 const int *hrow, const double *colels,
146 int nrows, int testvals,
147 const char *ROW, const char *COL)
148{
149 for (int irow=0; irow<nrows; irow++) {
150 if (hinrow[irow] > 0) {
151 CoinBigIndex krs = mrstrt[irow];
152 CoinBigIndex kre = krs + hinrow[irow];
153
154 for (CoinBigIndex k=krs; k<kre; k++) {
155 int jcol = hcol[k];
156 CoinBigIndex kcs = mcstrt[jcol];
157 CoinBigIndex kce = kcs + hincol[jcol];
158
159 CoinBigIndex kk = presolve_find_row1(irow, kcs, kce, hrow);
160 if (kk == kce) {
161 printf("MATRIX INCONSISTENT: can't find %s %d in %s %d\n",
162 ROW, irow, COL, jcol);
163 fflush(stdout);
164 abort();
165 }
166 if (testvals && colels[kk] != rowels[k]) {
167 printf("MATRIX INCONSISTENT: value for %s %d and %s %d\n",
168 ROW, irow, COL, jcol);
169 fflush(stdout);
170 abort();
171 }
172 }
173 }
174 }
175}
176#endif
177} // end unnamed file-local namespace
178
179/*
180 Utilizes matrix_consistent to check for equivalence of the column- and
181 row-major representations. Checks for presence of coefficients in the
182 column-major matrix, given presence in the row-major matrix, then checks
183 for presence in the row-major matrix given presence in the column-major
184 matrix. If testvals == true (default), the check also tests that the
185 coefficients have equal value.
186
187 See further comments with matrix_consistent.
188*/
189
190# if PRESOLVE_CONSISTENCY
191void presolve_consistent(const CoinPresolveMatrix *preObj, bool testvals)
192{
193 matrix_consistent(preObj->mrstrt_,preObj->hinrow_,preObj->hcol_,
194 preObj->rowels_,
195 preObj->mcstrt_,preObj->hincol_,preObj->hrow_,
196 preObj->colels_,
197 preObj->nrows_,testvals,"row","col") ;
198 matrix_consistent(preObj->mcstrt_,preObj->hincol_,preObj->hrow_,
199 preObj->colels_,
200 preObj->mrstrt_,preObj->hinrow_,preObj->hcol_,
201 preObj->rowels_,
202 preObj->ncols_,testvals,"col","row") ;
203}
204
205/*
206 Check the column- and/or row-major matrices for duplicates. By default, both
207 will be checked.
208*/
209
210void presolve_no_dups (const CoinPresolveMatrix *preObj,
211 bool doCol, bool doRow)
212
213{
214 if (doCol)
215 { no_majvec_dups(0,preObj->mcstrt_,preObj->hrow_,
216 preObj->hincol_,preObj->ncols_) ; }
217 if (doRow)
218 { no_majvec_dups(0,preObj->mrstrt_,preObj->hcol_,
219 preObj->hinrow_,preObj->nrows_) ; }
220
221 return ; }
222
223
224/*
225 As the name implies: scan for explicit zeros. By default, both matrices are
226 scanned.
227*/
228void presolve_no_zeros (const CoinPresolveMatrix *preObj,
229 bool doCol, bool doRow)
230{
231 if (doCol)
232 { check_majvec_nozeros(preObj->mcstrt_,preObj->colels_,preObj->hincol_,
233 preObj->ncols_) ; }
234 if (doRow)
235 { check_majvec_nozeros(preObj->mrstrt_,preObj->rowels_,preObj->hinrow_,
236 preObj->nrows_) ; }
237
238 return ; }
239
240/*
241 Lazy check on column lengths. Scan the row index array for the column.
242 If the relevant row length in the row-major rep is non-zero, assume we're ok.
243
244 Not advertised in CoinPresolvePsdebug.hpp.
245*/
246void presolve_hincol_ok(const int *mcstrt, const int *hincol,
247 const int *hinrow,
248 const int *hrow, int ncols)
249{
250 int jcol;
251
252 for (jcol=0; jcol<ncols; jcol++)
253 if (hincol[jcol] > 0) {
254 int kcs = mcstrt[jcol];
255 int kce = kcs + hincol[jcol];
256 int n=0;
257
258 int k;
259 for (k=kcs; k<kce; k++) {
260 int row = hrow[k];
261 if (hinrow[row] > 0)
262 n++;
263 }
264 if (n != hincol[jcol])
265 abort();
266 }
267}
268
269/*
270 Integrity checks for the linked lists that indicate major vector ordering
271 in the bulk storage area (minor index and coefficient arrays).
272 */
273void presolve_links_ok (const CoinPresolveMatrix *preObj,
274 bool doCol, bool doRow)
275{
276 if (doCol)
277 { links_ok(preObj->clink_,preObj->mcstrt_,
278 preObj->hincol_,preObj->ncols_) ; }
279 if (doRow)
280 { links_ok(preObj->rlink_,preObj->mrstrt_,
281 preObj->hinrow_,preObj->nrows_) ; }
282
283 return ; }
284
285
286
287/*
288 Routines to check a threaded matrix from a CoinPostsolve object.
289*/
290
291/*
292 Check that the column length agrees with the column thread. There must be
293 the correct number of coefficients, and the thread must end with the NO_LINK
294 marker.
295*/
296void presolve_check_threads (const CoinPostsolveMatrix *obj)
297
298{
299
300 CoinBigIndex *mcstrt = obj->mcstrt_ ;
301 int *hincol = obj->hincol_ ;
302 CoinBigIndex *link = obj->link_ ;
303 char *cdone = obj->cdone_ ;
304
305 int n = obj->ncols0_ ;
306
307/*
308 Scan the columns, checking only the ones that have been processed into the
309 constraint matrix.
310*/
311 for (int j = 0 ; j < n ; j++)
312 { if (!cdone[j]) continue ;
313
314 int lenj = hincol[j] ;
315 int k ;
316 for (k = mcstrt[j] ; k != NO_LINK && lenj > 0 ; k = link[k])
317 { assert(k >= 0 && k < obj->maxlink_) ;
318 lenj-- ; }
319
320 assert(k == NO_LINK && lenj == 0) ; }
321
322 return ; }
323
324/*
325 Check the free list. We're looking for gross corruption here. The notion is
326 that the free list plus elements in the matrix should add up to the capacity
327 of the bulk store.
328*/
329
330void presolve_check_free_list (const CoinPostsolveMatrix *obj, bool chkElemCnt)
331
332{
333
334 CoinBigIndex k = obj->free_list_ ;
335 CoinBigIndex freeCnt = 0 ;
336 CoinBigIndex maxlink = obj->maxlink_ ;
337 CoinBigIndex *link = obj->link_ ;
338/*
339 Redundancy in the data structure. These should always be equal.
340*/
341 assert(maxlink == obj->bulk0_) ;
342/*
343 Walk the free list portion of link. We should never point outside the bulk
344 store. If we ever come across an entry that's less than 0, it had better be
345 NO_LINK, the end marker.
346*/
347 while (k >= 0)
348 { assert(k < maxlink) ;
349 freeCnt++ ;
350 k = link[k] ; }
351 assert(k == NO_LINK) ;
352/*
353 And a final test: elements in the matrix plus free space should equal the
354 size of the bulk area. A good thought, but less than practical. Currently
355 postsolve doesn't track the number of elements in the matrix. But you might
356 find it useful if you're checking a newly constructed postsolve matrix. Even
357 then, you need to make sure nelems_ is correct. In the normal scheme of
358 things, this requires that somewhere there's a count of elements. Right now,
359 drop_empty_cols_action::presolve does this count, and you can get an accurate
360 value from the presolve object. assignPresolveToPostsolve will transfer this
361 value. Otherwise you're on your own --- your constructor must somehow find
362 this count. Using a standard CoinPackedMatrix is another way to get a count.
363*/
364 if (chkElemCnt)
365 { assert(obj->nelems_+freeCnt == maxlink) ; }
366
367
368 return ; }
369#endif
370
371
372/*
373 Routines to check solution and basis composition.
374*/
375
376/*
377 CoinPostsolveMatrix
378
379 This routine performs two checks on reduced costs held in rcosts_:
380 * The value held in rcosts_ is checked against the status of the
381 variable. Errors reported as "Bad rcost"
382 * The reduced cost is calculated from scratch and compared to the
383 value held in rcosts_. Errors reported as "Inacc rcost"
384
385 Remember that postsolve has a schizophrenic attitude about maximisation. All
386 transforms assume minimisation, and that's reflected in the reduced costs we
387 see here. And you must load duals and reduced costs with the correct sign for
388 minimisation. But, as a small courtesy (and a big inconsistency), postsolve
389 will negate objective coefficients for you. Hence the rather odd use of
390 maxmin.
391
392 The routine is specific to CoinPostsolveMatrix because the reduced cost
393 calculation requires traversal of (threaded) matrix columns.
394
395 NOTE: This routine holds static variables. It will detect when the problem
396 size changes and reinitialise. If you use presolve debugging over
397 multiple problems and you want to be dead sure of reinitialisation,
398 use the call presolve_check_reduced_costs(0), which will reinitialise
399 and return.
400*/
401# if PRESOLVE_DEBUG
402
403void presolve_check_reduced_costs (const CoinPostsolveMatrix *postObj)
404{
405
406 static bool warned = false ;
407 static double *warned_rcosts = 0 ;
408 static int allocSize = 0 ;
409 static const CoinPostsolveMatrix *lastObj = 0 ;
410
411/*
412 Is the client asking for reinitialisation only?
413*/
414 if (postObj == 0)
415 { warned = false ;
416 if (warned_rcosts != 0)
417 { delete[] warned_rcosts ;
418 warned_rcosts = 0 ; }
419 allocSize = 0 ;
420 lastObj = 0 ;
421 return ; }
422/*
423 *Should* the client have asked for reinitialisation?
424*/
425 int ncols0 = postObj->ncols0_ ;
426 if (allocSize < ncols0 || postObj != lastObj)
427 { warned = false ;
428 delete[] warned_rcosts ;
429 warned_rcosts = 0 ;
430 allocSize = 0 ;
431 lastObj = postObj ; }
432
433
434 double *rcosts = postObj->rcosts_ ;
435
436/*
437 By tracking values in warned_rcosts, we can produce a single message the
438 first time a value is determined to be incorrect.
439*/
440 if (!warned)
441 { warned = true ;
442 std::cout
443 << "reduced cost" << std::endl ;
444 warned_rcosts = new double[ncols0] ;
445 CoinZeroN(warned_rcosts,ncols0) ; }
446
447 double *colels = postObj->colels_ ;
448 int *hrow = postObj->hrow_ ;
449 int *mcstrt = postObj->mcstrt_ ;
450 int *hincol = postObj->hincol_ ;
451 CoinBigIndex *link = postObj->link_ ;
452
453 double *clo = postObj->clo_ ;
454 double *cup = postObj->cup_ ;
455
456 double *dcost = postObj->cost_ ;
457
458 double *sol = postObj->sol_ ;
459
460 char *cdone = postObj->cdone_ ;
461 char *rdone = postObj->rdone_ ;
462
463 const double ztoldj = postObj->ztoldj_ ;
464 const double ztolzb = postObj->ztolzb_ ;
465
466 double *rowduals = postObj->rowduals_ ;
467
468 double maxmin = postObj->maxmin_ ;
469 std::string strMaxmin((maxmin < 0)?"max":"min") ;
470 int checkCol=-1;
471/*
472 Scan all columns, but only check the ones that are marked as having been
473 postprocessed.
474*/
475 for (int j = 0 ; j < ncols0 ; j++)
476 { if (cdone[j] == 0) continue ;
477 const char *statjstr = postObj->columnStatusString(j) ;
478/*
479 Check the stored reduced cost for accuracy. See note above w.r.t. maxmin.
480*/
481 double dj = rcosts[j] ;
482 double wrndj = warned_rcosts[j] ;
483
484 { int ndx ;
485 CoinBigIndex k = mcstrt[j] ;
486 int len = hincol[j] ;
487 double chkdj = maxmin*dcost[j] ;
488 if (j==checkCol)
489 printf("dj for %d is %g - cost is %g\n",
490 j,dj,chkdj);
491 for (ndx = 0 ; ndx < len ; ndx++)
492 { int row = hrow[k] ;
493 PRESOLVEASSERT(rdone[row] != 0) ;
494 chkdj -= rowduals[row]*colels[k] ;
495 if (j==checkCol)
496 printf("row %d coeff %g dual %g => dj %g\n",
497 row,colels[k],rowduals[row],chkdj);
498
499 k = link[k] ; }
500 if (fabs(dj-chkdj) > ztoldj && wrndj != dj)
501 { std::cout
502 << "Inacc rcost: " << j << " " << statjstr << " "
503 << strMaxmin << " have " << dj
504 << " should be " << chkdj << " err " << fabs(dj-chkdj)
505 << std::endl ; } }
506/*
507 Check the stored reduced cost for consistency with the variable's status.
508 The cases are
509 * basic: (reduced cost) == 0
510 * at upper bound and not at lower bound: (reduced cost)*(maxmin) <= 0
511 * at lower bound and not at upper bound: (reduced cost)*(maxmin) >= 0
512 * not at either bound: any sign is correct (the variable can move either
513 way) but superbasic status is sufficiently exotic that it always
514 deserves a message. (There should be no superbasic variables at the
515 completion of postsolve.)
516 As a courtesy, show the reduced cost with the proper sign.
517*/
518 { double xj = sol[j] ;
519 double lj = clo[j] ;
520 double uj = cup[j] ;
521
522 if (postObj->columnIsBasic(j))
523 { if (fabs(dj) > ztoldj && wrndj != dj)
524 { std::cout
525 << "Bad rcost: " << j << " " << maxmin*dj
526 << " " << statjstr << " " << strMaxmin << std::endl ; } }
527 else
528 if (fabs(xj-uj) < ztolzb && fabs(xj-lj) > ztolzb)
529 { if (dj >= ztoldj && wrndj != dj)
530 { std::cout
531 << "Bad rcost: " << j << " " << maxmin*dj
532 << " " << statjstr << " " << strMaxmin << std::endl ; } }
533 else
534 if (fabs(xj-lj) < ztolzb && fabs(xj-uj) > ztolzb)
535 { if (dj <= -ztoldj && wrndj != dj)
536 { std::cout
537 << "Bad rcost: " << j << " " << maxmin*dj
538 << " " << statjstr << " " << strMaxmin << std::endl ; } }
539 else
540 if (fabs(xj-lj) > ztolzb && fabs(xj-uj) > ztolzb)
541 { if (fabs(dj) > ztoldj && wrndj != dj)
542 { std::cout
543 << "Superbasic rcost: " << j << " " << maxmin*dj
544 << " " << statjstr << " " << strMaxmin
545 << " lb "<< lj << " val " << xj << " ub "<< uj << std::endl ; } }
546 }
547
548 warned_rcosts[j] = rcosts[j] ; }
549
550}
551
552/*
553 CoinPostsolveMatrix
554
555 This routine checks the value and status of the dual variables. It
556 checks that the value and status of the dual agree with the row activity.
557 Errors are reported as "Bad dual"
558
559 See presolve_check_reduced_costs for an explanation of the use of maxmin.
560
561 Specific to CoinPostsolveMatrix due to the use of rdone. This could be fixed,
562 but probably better to clone the function and specialise it for
563 CoinPresolveMatrix.
564*/
565
566void presolve_check_duals (const CoinPostsolveMatrix *postObj)
567{
568
569
570 int nrows0 = postObj->nrows0_ ;
571
572 double *rowduals = postObj->rowduals_ ;
573
574 double *acts = postObj->acts_ ;
575 double *rup = postObj->rup_ ;
576 double *rlo = postObj->rlo_ ;
577
578 char *rdone = postObj->rdone_ ;
579
580 const double ztoldj = postObj->ztoldj_ ;
581 const double ztolzb = postObj->ztolzb_ ;
582
583 double maxmin = postObj->maxmin_ ;
584 std::string strMaxmin((maxmin < 0)?"max":"min") ;
585
586/*
587 Scan all processed rows. The rules are as for normal reduced costs, but
588 we need to remember the various flips and inversions. In summary, the correct
589 situation at optimality (minimisation) is:
590 * acts[i] == rup[i] ==> artificial NBLB ==> dual[i] < 0
591 * acts[i] == rlo[i] ==> artificial NBUB ==> dual[i] > 0
592
593 We can't say much about the dual for an equality. It can go either way. As a
594 courtesy, show the dual with the proper sign.
595*/
596 for (int i = 0 ; i < nrows0 ; i++)
597 { if (rdone[i] == 0) continue ;
598
599 double ui = rup[i] ;
600 double li = rlo[i] ;
601
602 if (ui-li < 1.0e-6) continue ;
603
604 double yi = rowduals[i] ;
605 double lhsi = acts[i] ;
606 const char *statistr = postObj->rowStatusString(i) ;
607
608
609 if (fabs(lhsi-li) < ztolzb)
610 { if (yi < -ztoldj)
611 { std::cout
612 << "Bad dual: " << i << " " << maxmin*yi
613 << " " << statistr << " " << strMaxmin << std::endl ; } }
614 else
615 if (fabs(lhsi-ui) < ztolzb)
616 { if (yi > ztoldj)
617 { std::cout
618 << "Bad dual: " << i << " " << maxmin*yi
619 << " " << statistr << " " << strMaxmin << std::endl ; } }
620 else
621 if (li < lhsi && lhsi < ui)
622 { if (fabs(yi) > ztoldj)
623 { std::cout
624 << "Bad dual: " << i << " " << maxmin*yi
625 << " " << statistr << " " << strMaxmin << std::endl ; } } }
626 return ; }
627
628
629
630/*
631 CoinPresolveMatrix
632
633 This routine will check the primal (column) solution for feasibility and
634 status. If there's no column solution (sol_), the routine bails out. If the
635 column solution is present, all else is assumed to be present.
636
637 chkColSol: check colum solution (primal variables)
638 0 - checks off
639 1 - check for NaN/Inf
640 *2 - check for above/below column bounds
641 chkRowAct: check row solution (evaluate constraint lhs)
642 0 - checks off
643 *1 - check for NaN/Inf
644 2 - check for inaccuracy, above/below row bounds
645 chkStatus: check for valid status of architectural variables
646 0 - checks off
647 *1 - checks on, if colstat_ exists
648
649 In general, the presolve transforms are not prepared to properly adjust the
650 row activity (reported as `Inacc RSOL'). Postsolve transforms do better. On
651 the bright side, the code seems to work just fine without maintaining row
652 activity. You probably don't want to use the level 2 checks for the row
653 solution, particularly in presolve.
654
655 To do: implement row status checks.
656
657 With a bit of thought, the various checks could be more cleanly separated
658 to require only the minimum information for each check.
659*/
660void presolve_check_sol (const CoinPresolveMatrix *preObj,
661 int chkColSol, int chkRowAct, int chkStatus)
662
663{
664 double *colels = preObj->colels_ ;
665 int *hrow = preObj->hrow_ ;
666 int *mcstrt = preObj->mcstrt_ ;
667 int *hincol = preObj->hincol_ ;
668 int *hinrow = preObj->hinrow_ ;
669
670 int n = preObj->ncols_ ;
671 int m = preObj->nrows_ ;
672
673/*
674 If there's no column solution, bail out now.
675*/
676 if (preObj->sol_ == 0) return ;
677
678 double *csol = preObj->sol_ ;
679 double *acts = preObj->acts_ ;
680 double *clo = preObj->clo_ ;
681 double *cup = preObj->cup_ ;
682 double *rlo = preObj->rlo_ ;
683 double *rup = preObj->rup_ ;
684
685 double tol = preObj->ztolzb_ ;
686
687 double *rsol = 0 ;
688 if (chkRowAct)
689 { rsol = new double[m] ;
690 memset(rsol,0,m*sizeof(double)) ; }
691
692/*
693 Open a loop to scan each column. For each column, do the following:
694 * Update the row solution (lhs value) by adding the contribution from
695 this column.
696 * Check for bogus values (NaN, infinity)
697 * Check for feasibility (value within column bounds)
698 * Check that the status of the variable agrees with the value and with the
699 lower and upper bounds. Free should have no bounds, superbasic should
700 have at least one.
701*/
702 for (int j = 0 ; j < n ; ++j)
703 { CoinBigIndex v = mcstrt[j] ;
704 int colLen = hincol[j] ;
705 double xj = csol[j] ;
706 double lj = clo[j] ;
707 double uj = cup[j] ;
708
709 if (chkRowAct)
710 { for (int u = 0 ; u < colLen ; ++u)
711 { int i = hrow[v] ;
712 double aij = colels[v] ;
713 v++ ;
714 rsol[i] += aij*xj ; } }
715
716 if (chkColSol&((1<<1)|(1<<0)))
717 { if (CoinIsnan(xj))
718 { printf("NaN CSOL: %d : lb = %g x = %g ub = %g\n",j,lj,xj,uj) ; }
719 if (xj <= -PRESOLVE_INF || xj >= PRESOLVE_INF)
720 { printf("Inf CSOL: %d : lb = %g x = %g ub = %g\n",j,lj,xj,uj) ; }
721 if (chkColSol > 1)
722 { if (xj < lj-tol)
723 { printf("low CSOL: %d : lb = %g x = %g ub = %g\n",j,lj,xj,uj) ; }
724 else
725 if (xj > uj+tol)
726 { printf("high CSOL: %d : lb = %g x = %g ub = %g\n",
727 j,lj,xj,uj) ; } } }
728 if (chkStatus && preObj->colstat_)
729 { CoinPrePostsolveMatrix::Status statj = preObj->getColumnStatus(j) ;
730 switch (statj)
731 { case CoinPrePostsolveMatrix::atUpperBound:
732 { if (uj >= PRESOLVE_INF || fabs(xj-uj) > tol)
733 { printf("Bad status CSOL: %d : status atUpperBound : ",j) ;
734 printf("lb = %g x = %g ub = %g\n",lj,xj,uj) ; }
735 break ; }
736 case CoinPrePostsolveMatrix::atLowerBound:
737 { if (lj <= -PRESOLVE_INF || fabs(xj-lj) > tol)
738 { printf("Bad status CSOL: %d : status atLowerBound : ",j) ;
739 printf("lb = %g x = %g ub = %g\n",lj,xj,uj) ; }
740 break ; }
741 case CoinPrePostsolveMatrix::isFree:
742 { if (lj > -PRESOLVE_INF || uj < PRESOLVE_INF)
743 { printf("Bad status CSOL: %d : status isFree : ",j) ;
744 printf("lb = %g x = %g ub = %g\n",lj,xj,uj) ; }
745 break ; }
746 case CoinPrePostsolveMatrix::superBasic:
747 { if (!(lj > -PRESOLVE_INF || uj < PRESOLVE_INF))
748 { printf("Bad status CSOL: %d : status superBasic : ",j) ;
749 printf("lb = %g x = %g ub = %g\n",lj,xj,uj) ; }
750 break ; }
751 case CoinPrePostsolveMatrix::basic:
752 { /* Nothing to do here. */
753 break ; }
754 default:
755 { printf("Bad status CSOL: %d : status unrecognized : ",j) ;
756 break ; } } } }
757/*
758 Now check the row solution. acts[i] is what presolve thinks we have, rsol[i]
759 is what we've just calculated while scanning the columns. We need only
760 check nontrivial rows (i.e., rows with length > 0). For each row,
761 * Check for bogus values (NaN, infinity)
762 * Check for accuracy (acts == rsol)
763 * Check for feasibility (rsol within row bounds)
764*/
765 tol *=1.0e3;
766 if (chkRowAct)
767 { for (int i = 0 ; i < m ; ++i)
768 { if (hinrow[i])
769 { double lhsi = acts[i] ;
770 double evali = rsol[i] ;
771 double li = rlo[i] ;
772 double ui = rup[i] ;
773
774 if (CoinIsnan(evali) || CoinIsnan(lhsi))
775 { printf("NaN RSOL: %d : lb = %g eval = %g (expected %g) ub = %g\n",
776 i,li,evali,lhsi,ui) ; }
777 if (evali <= -PRESOLVE_INF || evali >= PRESOLVE_INF ||
778 lhsi <= -PRESOLVE_INF || lhsi >= PRESOLVE_INF)
779 { printf("Inf RSOL: %d : lb = %g eval = %g (expected %g) ub = %g\n",
780 i,li,evali,lhsi,ui) ; }
781 if (chkRowAct > 1)
782 { if (fabs(evali-lhsi) > tol)
783 { printf("Inacc RSOL: %d : lb = %g eval = %g (expected %g) ub = %g\n",
784 i,li,evali,lhsi,ui) ; }
785 if (evali < li-tol || lhsi < li-tol)
786 { printf("low RSOL: %d : lb = %g eval = %g (expected %g) ub = %g\n",
787 i,li,evali,lhsi,ui) ; }
788 else
789 if (evali > ui+tol || lhsi > ui+tol)
790 { printf("high RSOL: %d : lb = %g eval = %g (expected %g) ub = %g\n",
791 i,li,evali,lhsi,ui) ; } } } }
792 delete [] rsol ; }
793 return ; }
794
795/*
796 CoinPostsolveMatrix
797
798 check_sol overload for CoinPostsolveMatrix. Parameters and functionality
799 identical to check_sol immediately above, but we have to remember we're
800 working with a threaded column-major representation.
801*/
802void presolve_check_sol (const CoinPostsolveMatrix *postObj,
803 int chkColSol, int chkRowAct, int chkStatus)
804
805{
806 double *colels = postObj->colels_ ;
807 int *hrow = postObj->hrow_ ;
808 int *mcstrt = postObj->mcstrt_ ;
809 int *hincol = postObj->hincol_ ;
810 int *link = postObj->link_ ;
811
812 int n = postObj->ncols_ ;
813 int m = postObj->nrows_ ;
814
815 double *csol = postObj->sol_ ;
816 double *acts = postObj->acts_ ;
817 double *clo = postObj->clo_ ;
818 double *cup = postObj->cup_ ;
819 double *rlo = postObj->rlo_ ;
820 double *rup = postObj->rup_ ;
821
822 double tol = postObj->ztolzb_ ;
823
824 double *rsol = 0 ;
825 if (chkRowAct)
826 { rsol = new double[m] ;
827 memset(rsol,0,m*sizeof(double)) ; }
828
829/*
830 Open a loop to scan each column. For each column, do the following:
831 * Update the row solution (lhs value) by adding the contribution from
832 this column.
833 * Check for bogus values (NaN, infinity)
834 * check that the status of the variable agrees with the value and with the
835 lower and upper bounds. Free should have no bounds, superbasic should
836 have at least one.
837*/
838 for (int j = 0 ; j < n ; ++j)
839 { CoinBigIndex v = mcstrt[j] ;
840 int colLen = hincol[j] ;
841 double xj = csol[j] ;
842 double lj = clo[j] ;
843 double uj = cup[j] ;
844
845 if (chkRowAct)
846 { for (int u = 0 ; u < colLen ; ++u)
847 { int i = hrow[v] ;
848 double aij = colels[v] ;
849 v = link[v] ;
850 rsol[i] += aij*xj ; } }
851 if (chkColSol)
852 { if (CoinIsnan(xj))
853 { printf("NaN CSOL: %d : lb = %g x = %g ub = %g\n",j,lj,xj,uj) ; }
854 if (xj <= -PRESOLVE_INF || xj >= PRESOLVE_INF)
855 { printf("Inf CSOL: %d : lb = %g x = %g ub = %g\n",j,lj,xj,uj) ; }
856 if (chkColSol > 1)
857 { if (xj < lj-tol)
858 { printf("low CSOL: %d : lb = %g x = %g ub = %g\n",j,lj,xj,uj) ; }
859 else
860 if (xj > uj+tol)
861 { printf("high CSOL: %d : lb = %g x = %g ub = %g\n",
862 j,lj,xj,uj) ; } } }
863 if (chkStatus && postObj->colstat_)
864 { CoinPrePostsolveMatrix::Status statj = postObj->getColumnStatus(j) ;
865 switch (statj)
866 { case CoinPrePostsolveMatrix::atUpperBound:
867 { if (uj >= PRESOLVE_INF || fabs(xj-uj) > tol)
868 { printf("Bad status CSOL: %d : status atUpperBound : ",j) ;
869 printf("lb = %g x = %g ub = %g\n",lj,xj,uj) ; }
870 break ; }
871 case CoinPrePostsolveMatrix::atLowerBound:
872 { if (lj <= -PRESOLVE_INF || fabs(xj-lj) > tol)
873 { printf("Bad status CSOL: %d : status atLowerBound : ",j) ;
874 printf("lb = %g x = %g ub = %g\n",lj,xj,uj) ; }
875 break ; }
876 case CoinPrePostsolveMatrix::isFree:
877 { if (lj > -PRESOLVE_INF || uj < PRESOLVE_INF)
878 { printf("Bad status CSOL: %d : status isFree : ",j) ;
879 printf("lb = %g x = %g ub = %g\n",lj,xj,uj) ; }
880 break ; }
881 case CoinPrePostsolveMatrix::superBasic:
882 { if (!(lj > -PRESOLVE_INF || uj < PRESOLVE_INF))
883 { printf("Bad status CSOL: %d : status superBasic : ",j) ;
884 printf("lb = %g x = %g ub = %g\n",lj,xj,uj) ; }
885 break ; }
886 case CoinPrePostsolveMatrix::basic:
887 { /* Nothing to do here. */
888 break ; }
889 default:
890 { printf("Bad status CSOL: %d : status unrecognized : ",j) ;
891 break ; } } } }
892/*
893 Now check the row solution. acts[i] is what presolve thinks we have, rsol[i]
894 is what we've just calculated while scanning the columns. CoinPostsolveMatrix
895 does not contain hinrow_, so we can't check for trivial rows (cf. check_sol
896 for CoinPresolveMatrix). For each row,
897 * Check for bogus values (NaN, infinity)
898*/
899 tol *= 1.0e4;
900 if (chkRowAct)
901 { for (int i = 0 ; i < m ; ++i)
902 { double lhsi = acts[i] ;
903 double evali = rsol[i] ;
904 double li = rlo[i] ;
905 double ui = rup[i] ;
906
907 if (CoinIsnan(evali) || CoinIsnan(lhsi))
908 { printf("NaN RSOL: %d : lb = %g eval = %g (expected %g) ub = %g\n",
909 i,li,evali,lhsi,ui) ; }
910 if (evali <= -PRESOLVE_INF || evali >= PRESOLVE_INF ||
911 lhsi <= -PRESOLVE_INF || lhsi >= PRESOLVE_INF)
912 { printf("Inf RSOL: %d : lb = %g eval = %g (expected %g) ub = %g\n",
913 i,li,evali,lhsi,ui) ; }
914 if (chkRowAct > 1)
915 { if (fabs(evali-lhsi) > tol)
916 { printf("Inacc RSOL: %d : lb = %g eval = %g (expected %g) ub = %g\n",
917 i,li,evali,lhsi,ui) ; }
918 if (evali < li-tol || lhsi < li-tol)
919 { printf("low RSOL: %d : lb = %g eval = %g (expected %g) ub = %g\n",
920 i,li,evali,lhsi,ui) ; }
921 else
922 if (evali > ui+tol || lhsi > ui+tol)
923 { printf("high RSOL: %d : lb = %g eval = %g (expected %g) ub = %g\n",
924 i,li,evali,lhsi,ui) ; } } }
925 delete [] rsol ; }
926 return ; }
927
928/*
929 CoinPostsolveMatrix
930
931 Make sure that the number of basic variables is correct.
932*/
933void presolve_check_nbasic (const CoinPostsolveMatrix *postObj)
934
935{
936
937 int ncols0 = postObj->ncols0_ ;
938 int nrows0 = postObj->nrows0_ ;
939
940 char *cdone = postObj->cdone_ ;
941 char *rdone = postObj->rdone_ ;
942
943 int nbasic = 0 ;
944 int ncdone = 0;
945 int nrdone = 0;
946 int ncb = 0;
947 int nrb = 0;
948
949 for (int j = 0 ; j < ncols0 ; j++)
950 {
951 if (cdone[j] != 0 && postObj->columnIsBasic(j))
952 { nbasic++ ;
953 ncb++ ; }
954 if (cdone[j])
955 ncdone++ ;
956 }
957
958 for (int i = 0 ; i < nrows0 ; i++)
959 {
960 if (rdone[i] && postObj->rowIsBasic(i))
961 { nbasic++ ;
962 nrb++ ; }
963 if (rdone[i])
964 nrdone++ ;
965 }
966
967 if (nbasic != postObj->nrows_)
968 { printf("WRONG NUMBER NBASIC: is %d, should be %d; ",
969 nbasic,postObj->nrows_) ;
970 printf("cdone %d, col basic %d, rdone %d, row basic %d.\n",
971 ncdone,ncb,nrdone,nrb) ;
972 fflush(stdout) ; }
973 return ; }
974
975
976/*
977 CoinPresolveMatrix
978
979 Overload of presolve_check_nbasic for a CoinPresolveMatrix. There may not be
980 a solution, eh?
981*/
982void presolve_check_nbasic (const CoinPresolveMatrix *preObj)
983
984{
985
986 if (preObj->sol_ == 0) return ;
987
988 int ncols = preObj->ncols_ ;
989 int nrows = preObj->nrows_ ;
990
991 int nbasic = 0 ;
992 int ncb = 0;
993 int nrb = 0;
994
995 for (int j = 0 ; j < ncols ; j++)
996 {
997 if (preObj->columnIsBasic(j))
998 { nbasic++ ;
999 ncb++ ; }
1000 }
1001
1002 for (int i = 0 ; i < nrows ; i++)
1003 {
1004 if (preObj->rowIsBasic(i))
1005 { nbasic++ ;
1006 nrb++ ; }
1007 }
1008
1009 if (nbasic != nrows)
1010 { printf("WRONG NUMBER NBASIC: is: %d should be: %d;",
1011 nbasic,nrows) ;
1012 printf(" cb %d, rb %d.\n",ncb,nrb);
1013 fflush(stdout) ; }
1014 return ; }
1015
1016#endif
1017/*
1018 Original comment: I've forgotton what this is all about
1019
1020 Looks to me like it's confirming that the columns flagged as basic indeed
1021 have enough coefficients between them to cover the basis. It'd be serious
1022 work to get this going again. Waaaaaay out of date. -- lh, 040831 --
1023*/
1024# if 0
1025void check_pivots (const int *mrstrt, const int *hinrow, const int *hcol,
1026 int nrows, const unsigned char *colstat,
1027 const unsigned char *rowstat, int ncols)
1028{
1029 int i ;
1030 int nbasic = 0 ;
1031 int gotone = 1 ;
1032 int stillmore ;
1033
1034 return ;
1035
1036 int *bcol = new int[nrows] ;
1037 memset(bcol, -1, nrows*sizeof(int)) ;
1038
1039 char *coldone = new char[ncols] ;
1040 memset(coldone, 0, ncols) ;
1041
1042 while (gotone) {
1043 gotone = 0 ;
1044 stillmore = 0 ;
1045 for (i=0; i<nrows; i++)
1046 if (!postObj->rowIsBasic(i)) {
1047 int krs = mrstrt[i] ;
1048 int kre = mrstrt[i] + hinrow[i] ;
1049 int nb = 0 ;
1050 int kk ;
1051 for (int k=krs; k<kre; k++)
1052 if (postObj->columnIsBasic(hcol[k]) && !coldone[hcol[k]]) {
1053 nb++ ;
1054 kk = k ;
1055 if (nb > 1)
1056 break ;
1057 }
1058 if (nb == 1) {
1059 PRESOLVEASSERT(bcol[i] == -1) ;
1060 bcol[i] = hcol[kk] ;
1061 coldone[hcol[kk]] = 1 ;
1062 nbasic++ ;
1063 gotone = 1 ;
1064 }
1065 else
1066 stillmore = 1 ;
1067 }
1068 }
1069 PRESOLVEASSERT(!stillmore) ;
1070
1071 for (i=0; i<nrows; i++)
1072 if (postObj->rowIsBasic(i)) {
1073 int krs = mrstrt[i] ;
1074 int kre = mrstrt[i] + hinrow[i] ;
1075 for (int k=krs; k<kre; k++)
1076 PRESOLVEASSERT(!postObj->columnIsBasic(hcol[k]) || coldone[hcol[k]]) ;
1077 nbasic++ ;
1078 }
1079 PRESOLVEASSERT(nbasic == nrows) ;
1080}
1081
1082# endif
1083