source: trunk/PresolvePsdebug.cpp @ 55

Last change on this file since 55 was 55, checked in by forrest, 17 years ago

Making it a bit more likely Visual Studio will work

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.3 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include <new>
5#include <stdio.h>
6#include <math.h>
7
8#include "PresolveMatrix.hpp"
9
10
11static inline const double *
12lookup_by_col(const int *mcstrt, const int *hrow,
13              const int *hincol, const int *hinrow, const double *colels,
14              int row, int col)
15{
16  int kcs = mcstrt[col];
17  int kce = kcs + hincol[col];
18  int k;
19
20  if (hinrow[row] <= 0)
21    return (0);
22
23  for (k=kcs; k<kce; k++) {
24    int r = hrow[k];
25    if (r == row)
26      return (&colels[k]);
27  }
28  return (0);
29}
30
31static inline void
32no_dups(const char *done,
33        const int *mcstrt, const int *hrow, const int *hincol, int ncols)
34{
35#if     DEBUG_PRESOLVE
36  for (int jcol=0; jcol<ncols; jcol++) 
37    if ((!done || done[jcol]) && hincol[jcol] > 0) {
38      int kcs = mcstrt[jcol];
39      int kce = kcs + hincol[jcol];
40
41      for (int k=kcs; k<kce; k++) {
42        int row = hrow[k];
43
44        PRESOLVEASSERT(presolve_find_row1(row, k+1, kce, hrow) == kce);
45      }
46    }
47#endif
48}
49
50void presolve_no_zeros(const int *mcstrt, const double *colels, const int *hincol, int ncols)
51{
52#if     DEBUG_PRESOLVE
53  for (int jcol=0; jcol<ncols; jcol++) 
54    if (hincol[jcol] > 0) {
55      int kcs = mcstrt[jcol];
56      int kce = kcs + hincol[jcol];
57
58      for (int k=kcs; k<kce; k++) {
59        PRESOLVEASSERT(fabs(colels[k]) > ZTOLDP);
60      }
61    }
62#endif
63}
64
65
66void presolve_hincol_ok(const int *mcstrt, const int *hincol,
67               const int *hinrow,
68               const int *hrow, int ncols)
69{
70#if     CHECK_CONSISTENCY
71  int jcol;
72
73  for (jcol=0; jcol<ncols; jcol++) 
74    if (hincol[jcol] > 0) {
75      int kcs = mcstrt[jcol];
76      int kce = kcs + hincol[jcol];
77      int n=0;
78     
79      int k;
80      for (k=kcs; k<kce; k++) {
81        int row = hrow[k];
82        if (hinrow[row] > 0)
83          n++;
84      }
85      if (n != hincol[jcol])
86        abort();
87    }
88#endif
89}
90
91/*
92 * The linked list
93 */
94void presolve_links_ok(presolvehlink *link, int *starts, int *lengths, int n)
95{
96#if     CHECK_CONSISTENCY
97  int i;
98
99  for (i=0; i<n; i++) {
100    int pre = link[i].pre;
101    int suc = link[i].suc;
102
103    if (pre != NO_LINK) {
104      PRESOLVEASSERT(0 <= pre && pre <= n);
105      PRESOLVEASSERT(link[pre].suc == i);
106    }
107    if (suc != NO_LINK) {
108      PRESOLVEASSERT(0 <= suc && suc <= n);
109      PRESOLVEASSERT(link[suc].pre == i);
110    }
111  }
112
113  for (i=0; i<n; i++) 
114    if (link[i].pre == NO_LINK)
115      break;
116  PRESOLVEASSERT(i<n);
117
118  while (i != NO_LINK) {
119    if (link[i].suc != NO_LINK) 
120      PRESOLVEASSERT(starts[i] + lengths[i] <= starts[link[i].suc]);
121    i = link[i].suc;
122  }
123#endif
124}
125
126// I've forgotton what this is all about
127void check_pivots(const int *mrstrt, const int *hinrow, const int *hcol, int nrows,
128                  const unsigned char *colstat, const unsigned char *rowstat,
129                  int ncols)
130{
131#if 0
132  int i;
133  int nbasic = 0;
134  int gotone = 1;
135  int stillmore;
136
137  return;
138
139  int *bcol = new int[nrows];
140  memset(bcol, -1, nrows*sizeof(int));
141
142  char *coldone = new char[ncols];
143  memset(coldone, 0, ncols);
144
145  while (gotone) {
146    gotone = 0;
147    stillmore = 0;
148    for (i=0; i<nrows; i++)
149      if (!prob->rowIsBasic(i)) {
150        int krs = mrstrt[i];
151        int kre = mrstrt[i] + hinrow[i];
152        int nb = 0;
153        int kk;
154        for (int k=krs; k<kre; k++)
155          if (prob->columnIsBasic(hcol[k]) && !coldone[hcol[k]]) {
156            nb++;
157            kk = k;
158            if (nb > 1)
159              break;
160          }
161        if (nb == 1) {
162          PRESOLVEASSERT(bcol[i] == -1);
163          bcol[i] = hcol[kk];
164          coldone[hcol[kk]] = 1;
165          nbasic++;
166          gotone = 1;
167        }
168        else
169          stillmore = 1;
170      }
171  }
172  PRESOLVEASSERT(!stillmore);
173
174  for (i=0; i<nrows; i++)
175    if (prob->rowIsBasic(i)) {
176      int krs = mrstrt[i];
177      int kre = mrstrt[i] + hinrow[i];
178      for (int k=krs; k<kre; k++)
179        PRESOLVEASSERT(!prob->columnIsBasic(hcol[k]) || coldone[hcol[k]]);
180      nbasic++;
181    }
182  PRESOLVEASSERT(nbasic == nrows);
183#endif
184}
185
186
187
188static inline void a_ok(PostsolveMatrix *prob)
189{
190#if 0
191  static int warned = 0;
192  if (!warned) {
193    warned = 1;
194    printf("********** NO CHECKING!!!!!!! \n");
195  }
196
197  double *colels        = prob->colels;
198  int *hrow             = prob->hrow;
199  int *mcstrt           = prob->mcstrt;
200  int *hincol           = prob->hincol;
201
202  int *colstat  = prob->colstat;
203  int *rowstat  = prob->rowstat;
204
205  double *clo   = prob->clo;
206  double *cup   = prob->cup;
207
208  double *dcost = prob->cost;
209
210  double *sol   = prob->sol;
211  double *rcosts        = prob->rcosts;
212
213  char *cdone   = prob->cdone;
214  char *rdone   = prob->rdone;
215
216  const double ztoldj   = prob->ztoldj;
217  const double ztolzb   = prob->ztolzb;
218
219  double *rowduals = prob->rowduals;
220
221  int ncols0            = prob->ncols0;
222
223#if 0
224  {
225    int ncols           = prob->ncols;
226    int nrows           = prob->nrows;
227    int *mrstrt         = prob->mrstrt;
228    int *hinrow         = prob->hinrow;
229    int *hcol           = prob->hcol;
230
231    no_dups(cdone, mcstrt, hrow, hincol, ncols);
232    no_dups(rdone, mrstrt, hcol, hinrow, nrows);
233  }
234#endif
235
236  /* DEBUG HACK */
237  static double *warned_rcosts;
238  if (!warned_rcosts) {
239    warned_rcosts = new double[ncols0];
240
241    for (int i=0; i<ncols0; i++)
242      warned_rcosts[i] = rcosts[i];
243  }
244
245  for (int j=0; j<ncols0; j++)
246    if (cdone[j]) {
247      //int i = &paction[npaction-1]-pa;
248      int i = -666;
249
250      if (prob->columnIsBasic(j)) {
251        if (! (fabs(rcosts[j]) < ztoldj) &&
252            warned_rcosts[j] != rcosts[j])
253          printf("ITER[%d]: basic rcosts[%d]:  %g\n", i, j, rcosts[j]);
254      } else if (fabs(sol[j] - cup[j]) < ztolzb &&
255                 ! (fabs(sol[j] - clo[j]) < ztolzb)) {
256        if (! (rcosts[j] <= ztoldj) &&
257            warned_rcosts[j] != rcosts[j])
258          printf("ITER[%d]: ub rcosts[%d]:  %g\n", i, j, rcosts[j]);
259      } else if (fabs(sol[j] - clo[j]) < ztolzb &&
260                 ! (fabs(sol[j] - cup[j]) < ztolzb)){
261        if (! (rcosts[j] >= -ztoldj) &&
262            warned_rcosts[j] != rcosts[j])
263          printf("ITER[%d]: lb rcosts[%d]:  %g\n", i, j, rcosts[j]);
264      } else if (! (fabs(sol[j] - clo[j]) < ztolzb) &&
265                 ! (fabs(sol[j] - cup[j]) < ztolzb)) {
266        printf("SUPERBASIC (cost=%g):  %d %g < %g < %g\n",
267               dcost[j], j, clo[j], sol[j], cup[j]);
268      }
269
270      {
271        int kcs = mcstrt[j];
272        int kce = mcstrt[j] + hincol[j];
273        int k;
274        double dj = dcost[j];
275        int row0 = (kcs<kce ? hrow[kcs] : 0);
276
277        for (k=kcs; k<kce; k++) {
278          int row = hrow[k];
279          PRESOLVEASSERT(rdone[row]);
280          PRESOLVEASSERT(row != row0 || k==kcs);
281          if (rdone[row])
282            dj -= rowduals[row] * colels[k];
283        }
284        if (! (fabs(rcosts[j] - dj) < ztoldj) &&
285            warned_rcosts[j] != rcosts[j])
286          printf("ITER[%d]: rcosts[%d]:  %g <-> %g\n", i, j, rcosts[j], dj);
287      }
288    }
289
290  {
291    int i;
292    for (i=0; i<ncols0; i++)
293      warned_rcosts[i] = rcosts[i];
294  }
295#endif
296}
297
298
Note: See TracBrowser for help on using the repository browser.