source: branches/devel-1/PresolvePsdebug.cpp @ 117

Last change on this file since 117 was 49, checked in by ladanyi, 17 years ago

everything compiles

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