00001
00002
00003
00004 #include <new>
00005 #include <stdio.h>
00006 #include <math.h>
00007
00008 #include "CoinPresolveMatrix.hpp"
00009
00010
00011 static inline const double *
00012 lookup_by_col(const int *mcstrt, const int *hrow,
00013 const int *hincol, const int *hinrow, const double *colels,
00014 int row, int col)
00015 {
00016 int kcs = mcstrt[col];
00017 int kce = kcs + hincol[col];
00018 int k;
00019
00020 if (hinrow[row] <= 0)
00021 return (0);
00022
00023 for (k=kcs; k<kce; k++) {
00024 int r = hrow[k];
00025 if (r == row)
00026 return (&colels[k]);
00027 }
00028 return (0);
00029 }
00030
00031 static inline void
00032 no_dups(const char *done,
00033 const int *mcstrt, const int *hrow, const int *hincol, int ncols)
00034 {
00035 #if DEBUG_PRESOLVE
00036 for (int jcol=0; jcol<ncols; jcol++)
00037 if ((!done || done[jcol]) && hincol[jcol] > 0) {
00038 int kcs = mcstrt[jcol];
00039 int kce = kcs + hincol[jcol];
00040
00041 for (int k=kcs; k<kce; k++) {
00042 int row = hrow[k];
00043
00044 PRESOLVEASSERT(presolve_find_row1(row, k+1, kce, hrow) == kce);
00045 }
00046 }
00047 #endif
00048 }
00049
00050 void presolve_no_zeros(const int *mcstrt, const double *colels, const int *hincol, int ncols)
00051 {
00052 #if DEBUG_PRESOLVE
00053 for (int jcol=0; jcol<ncols; jcol++)
00054 if (hincol[jcol] > 0) {
00055 int kcs = mcstrt[jcol];
00056 int kce = kcs + hincol[jcol];
00057
00058 for (int k=kcs; k<kce; k++) {
00059 PRESOLVEASSERT(fabs(colels[k]) > ZTOLDP);
00060 }
00061 }
00062 #endif
00063 }
00064
00065
00066 void presolve_hincol_ok(const int *mcstrt, const int *hincol,
00067 const int *hinrow,
00068 const int *hrow, int ncols)
00069 {
00070 #if CHECK_CONSISTENCY
00071 int jcol;
00072
00073 for (jcol=0; jcol<ncols; jcol++)
00074 if (hincol[jcol] > 0) {
00075 int kcs = mcstrt[jcol];
00076 int kce = kcs + hincol[jcol];
00077 int n=0;
00078
00079 int k;
00080 for (k=kcs; k<kce; k++) {
00081 int row = hrow[k];
00082 if (hinrow[row] > 0)
00083 n++;
00084 }
00085 if (n != hincol[jcol])
00086 abort();
00087 }
00088 #endif
00089 }
00090
00091
00092
00093
00094 void presolve_links_ok(presolvehlink *link, int *starts, int *lengths, int n)
00095 {
00096 #if CHECK_CONSISTENCY
00097 int i;
00098
00099 for (i=0; i<n; i++) {
00100 int pre = link[i].pre;
00101 int suc = link[i].suc;
00102
00103 if (pre != NO_LINK) {
00104 PRESOLVEASSERT(0 <= pre && pre <= n);
00105 PRESOLVEASSERT(link[pre].suc == i);
00106 }
00107 if (suc != NO_LINK) {
00108 PRESOLVEASSERT(0 <= suc && suc <= n);
00109 PRESOLVEASSERT(link[suc].pre == i);
00110 }
00111 }
00112
00113 for (i=0; i<n; i++)
00114 if (link[i].pre == NO_LINK)
00115 break;
00116 PRESOLVEASSERT(i<n);
00117
00118 while (i != NO_LINK) {
00119 if (link[i].suc != NO_LINK)
00120 PRESOLVEASSERT(starts[i] + lengths[i] <= starts[link[i].suc]);
00121 i = link[i].suc;
00122 }
00123 #endif
00124 }
00125
00126
00127 void check_pivots(const int *mrstrt, const int *hinrow, const int *hcol, int nrows,
00128 const unsigned char *colstat, const unsigned char *rowstat,
00129 int ncols)
00130 {
00131 #if 0
00132 int i;
00133 int nbasic = 0;
00134 int gotone = 1;
00135 int stillmore;
00136
00137 return;
00138
00139 int *bcol = new int[nrows];
00140 memset(bcol, -1, nrows*sizeof(int));
00141
00142 char *coldone = new char[ncols];
00143 memset(coldone, 0, ncols);
00144
00145 while (gotone) {
00146 gotone = 0;
00147 stillmore = 0;
00148 for (i=0; i<nrows; i++)
00149 if (!prob->rowIsBasic(i)) {
00150 int krs = mrstrt[i];
00151 int kre = mrstrt[i] + hinrow[i];
00152 int nb = 0;
00153 int kk;
00154 for (int k=krs; k<kre; k++)
00155 if (prob->columnIsBasic(hcol[k]) && !coldone[hcol[k]]) {
00156 nb++;
00157 kk = k;
00158 if (nb > 1)
00159 break;
00160 }
00161 if (nb == 1) {
00162 PRESOLVEASSERT(bcol[i] == -1);
00163 bcol[i] = hcol[kk];
00164 coldone[hcol[kk]] = 1;
00165 nbasic++;
00166 gotone = 1;
00167 }
00168 else
00169 stillmore = 1;
00170 }
00171 }
00172 PRESOLVEASSERT(!stillmore);
00173
00174 for (i=0; i<nrows; i++)
00175 if (prob->rowIsBasic(i)) {
00176 int krs = mrstrt[i];
00177 int kre = mrstrt[i] + hinrow[i];
00178 for (int k=krs; k<kre; k++)
00179 PRESOLVEASSERT(!prob->columnIsBasic(hcol[k]) || coldone[hcol[k]]);
00180 nbasic++;
00181 }
00182 PRESOLVEASSERT(nbasic == nrows);
00183 #endif
00184 }
00185
00186
00187
00188 static inline void a_ok(CoinPostsolveMatrix *prob)
00189 {
00190 #if 0
00191 static int warned = 0;
00192 if (!warned) {
00193 warned = 1;
00194 printf("********** NO CHECKING!!!!!!! \n");
00195 }
00196
00197 double *colels = prob->colels;
00198 int *hrow = prob->hrow;
00199 int *mcstrt = prob->mcstrt;
00200 int *hincol = prob->hincol;
00201
00202 int *colstat = prob->colstat;
00203 int *rowstat = prob->rowstat;
00204
00205 double *clo = prob->clo;
00206 double *cup = prob->cup;
00207
00208 double *dcost = prob->cost;
00209
00210 double *sol = prob->sol;
00211 double *rcosts = prob->rcosts;
00212
00213 char *cdone = prob->cdone;
00214 char *rdone = prob->rdone;
00215
00216 const double ztoldj = prob->ztoldj;
00217 const double ztolzb = prob->ztolzb;
00218
00219 double *rowduals = prob->rowduals;
00220
00221 int ncols0 = prob->ncols0;
00222
00223 #if 0
00224 {
00225 int ncols = prob->ncols;
00226 int nrows = prob->nrows;
00227 int *mrstrt = prob->mrstrt;
00228 int *hinrow = prob->hinrow;
00229 int *hcol = prob->hcol;
00230
00231 no_dups(cdone, mcstrt, hrow, hincol, ncols);
00232 no_dups(rdone, mrstrt, hcol, hinrow, nrows);
00233 }
00234 #endif
00235
00236
00237 static double *warned_rcosts;
00238 if (!warned_rcosts) {
00239 warned_rcosts = new double[ncols0];
00240
00241 for (int i=0; i<ncols0; i++)
00242 warned_rcosts[i] = rcosts[i];
00243 }
00244
00245 for (int j=0; j<ncols0; j++)
00246 if (cdone[j]) {
00247
00248 int i = -666;
00249
00250 if (prob->columnIsBasic(j)) {
00251 if (! (fabs(rcosts[j]) < ztoldj) &&
00252 warned_rcosts[j] != rcosts[j])
00253 printf("ITER[%d]: basic rcosts[%d]: %g\n", i, j, rcosts[j]);
00254 } else if (fabs(sol[j] - cup[j]) < ztolzb &&
00255 ! (fabs(sol[j] - clo[j]) < ztolzb)) {
00256 if (! (rcosts[j] <= ztoldj) &&
00257 warned_rcosts[j] != rcosts[j])
00258 printf("ITER[%d]: ub rcosts[%d]: %g\n", i, j, rcosts[j]);
00259 } else if (fabs(sol[j] - clo[j]) < ztolzb &&
00260 ! (fabs(sol[j] - cup[j]) < ztolzb)){
00261 if (! (rcosts[j] >= -ztoldj) &&
00262 warned_rcosts[j] != rcosts[j])
00263 printf("ITER[%d]: lb rcosts[%d]: %g\n", i, j, rcosts[j]);
00264 } else if (! (fabs(sol[j] - clo[j]) < ztolzb) &&
00265 ! (fabs(sol[j] - cup[j]) < ztolzb)) {
00266 printf("SUPERBASIC (cost=%g): %d %g < %g < %g\n",
00267 dcost[j], j, clo[j], sol[j], cup[j]);
00268 }
00269
00270 {
00271 int kcs = mcstrt[j];
00272 int kce = mcstrt[j] + hincol[j];
00273 int k;
00274 double dj = dcost[j];
00275 int row0 = (kcs<kce ? hrow[kcs] : 0);
00276
00277 for (k=kcs; k<kce; k++) {
00278 int row = hrow[k];
00279 PRESOLVEASSERT(rdone[row]);
00280 PRESOLVEASSERT(row != row0 || k==kcs);
00281 if (rdone[row])
00282 dj -= rowduals[row] * colels[k];
00283 }
00284 if (! (fabs(rcosts[j] - dj) < ztoldj) &&
00285 warned_rcosts[j] != rcosts[j])
00286 printf("ITER[%d]: rcosts[%d]: %g <-> %g\n", i, j, rcosts[j], dj);
00287 }
00288 }
00289
00290 {
00291 int i;
00292 for (i=0; i<ncols0; i++)
00293 warned_rcosts[i] = rcosts[i];
00294 }
00295 #endif
00296 }
00297
00298