source: trunk/Clp/src/IdiSolve.cpp @ 1525

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 53.8 KB
Line 
1/* $Id: IdiSolve.cpp 1525 2010-02-26 17:27:59Z mjs $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4
5#include "CoinPragma.hpp"
6#include <stdio.h>
7#include <stdarg.h>
8#include <stdlib.h>
9#include <math.h>
10#include "CoinHelperFunctions.hpp"
11#include "Idiot.hpp"
12#define FIT
13#ifdef FIT
14#define HISTORY 8
15#else
16#define HISTORY 7
17#endif
18#define NSOLVE HISTORY-1
19static void solveSmall(int nsolve, double **aIn, double **a, double * b)
20{
21     int i, j;
22     /* copy */
23     for (i = 0; i < nsolve; i++) {
24          for (j = 0; j < nsolve; j++) {
25               a[i][j] = aIn[i][j];
26          }
27     }
28     for (i = 0; i < nsolve; i++) {
29          /* update using all previous */
30          double diagonal;
31          int j;
32          for (j = i; j < nsolve; j++) {
33               int k;
34               double value = a[i][j];
35               for (k = 0; k < i; k++) {
36                    value -= a[k][i] * a[k][j];
37               }
38               a[i][j] = value;
39          }
40          diagonal = a[i][i];
41          if (diagonal < 1.0e-20) {
42               diagonal = 0.0;
43          } else {
44               diagonal = 1.0 / sqrt(diagonal);
45          }
46          a[i][i] = diagonal;
47          for (j = i + 1; j < nsolve; j++) {
48               a[i][j] *= diagonal;
49          }
50     }
51     /* update */
52     for (i = 0; i < nsolve; i++) {
53          int j;
54          double value = b[i];
55          for (j = 0; j < i; j++) {
56               value -= b[j] * a[j][i];
57          }
58          value *= a[i][i];
59          b[i] = value;
60     }
61     for (i = nsolve - 1; i >= 0; i--) {
62          int j;
63          double value = b[i];
64          for (j = i + 1; j < nsolve; j++) {
65               value -= b[j] * a[i][j];
66          }
67          value *= a[i][i];
68          b[i] = value;
69     }
70}
71IdiotResult
72Idiot::objval(int nrows, int ncols, double * rowsol , double * colsol,
73              double * pi, double * /*djs*/, const double * cost ,
74              const double * /*rowlower*/,
75              const double * rowupper, const double * /*lower*/,
76              const double * /*upper*/, const double * elemnt,
77              const int * row, const CoinBigIndex * columnStart,
78              const int * length, int extraBlock, int * rowExtra,
79              double * solExtra, double * elemExtra, double * /*upperExtra*/,
80              double * costExtra, double weight)
81{
82     IdiotResult result;
83     double objvalue = 0.0;
84     double sum1 = 0.0, sum2 = 0.0;
85     int i;
86     for (i = 0; i < nrows; i++) {
87          rowsol[i] = -rowupper[i];
88     }
89     for (i = 0; i < ncols; i++) {
90          CoinBigIndex j;
91          double value = colsol[i];
92          if (value) {
93               objvalue += value * cost[i];
94               if (elemnt) {
95                    for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
96                         int irow = row[j];
97                         rowsol[irow] += elemnt[j] * value;
98                    }
99               } else {
100                    for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
101                         int irow = row[j];
102                         rowsol[irow] += value;
103                    }
104               }
105          }
106     }
107     /* adjust to make as feasible as possible */
108     /* no */
109     if (extraBlock) {
110          for (i = 0; i < extraBlock; i++) {
111               double element = elemExtra[i];
112               int irow = rowExtra[i];
113               objvalue += solExtra[i] * costExtra[i];
114               rowsol[irow] += solExtra[i] * element;
115          }
116     }
117     for (i = 0; i < nrows; i++) {
118          double value = rowsol[i];
119          sum1 += fabs(value);
120          sum2 += value * value;
121          pi[i] = -2.0 * weight * value;
122     }
123     result.infeas = sum1;
124     result.objval = objvalue;
125     result.weighted = objvalue + weight * sum2;
126     result.sumSquared = sum2;
127     return result;
128}
129IdiotResult
130Idiot::IdiSolve(
131     int nrows, int ncols, double * rowsol , double * colsol,
132     double * pi, double * djs, const double * origcost , double * rowlower,
133     double * rowupper, const double * lower,
134     const double * upper, const double * elemnt,
135     const int * row, const CoinBigIndex * columnStart,
136     const int * length, double * lambda,
137     int maxIts, double mu, double drop,
138     double maxmin, double offset,
139     int strategy, double djTol, double djExit, double djFlag,
140     CoinThreadRandom * randomNumberGenerator)
141{
142     IdiotResult result;
143     int  i, j, k, iter;
144     double value = 0.0, objvalue = 0.0, weightedObj = 0.0;
145     double tolerance = 1.0e-8;
146     double * history[HISTORY+1];
147     int ncolx;
148     int nChange;
149     int extraBlock = 0;
150     int * rowExtra = NULL;
151     double * solExtra = NULL;
152     double * elemExtra = NULL;
153     double * upperExtra = NULL;
154     double * costExtra = NULL;
155     double * useCostExtra = NULL;
156     double * saveExtra = NULL;
157     double * cost = NULL;
158     double saveValue = 1.0e30;
159     double saveOffset = offset;
160     double useOffset = offset;
161     /*#define NULLVECTOR*/
162#ifndef NULLVECTOR
163     int nsolve = NSOLVE;
164#else
165     int nsolve = NSOLVE + 1; /* allow for null vector */
166#endif
167     int nflagged;
168     double * thetaX;
169     double * djX;
170     double * bX;
171     double * vX;
172     double ** aX;
173     double **aworkX;
174     double ** allsum;
175     double * saveSol = 0;
176     const double * useCost = cost;
177     double bestSol = 1.0e60;
178     double weight = 0.5 / mu;
179     char * statusSave = new char[2*ncols];
180     char * statusWork = statusSave + ncols;
181#define DJTEST 5
182     double djSave[DJTEST];
183     double largestDj = 0.0;
184     double smallestDj = 1.0e60;
185     double maxDj = 0.0;
186     int doFull = 0;
187#define SAVEHISTORY 10
188#define EVERY (2*SAVEHISTORY)
189#define AFTER SAVEHISTORY*(HISTORY+1)
190#define DROP 5
191     double after = AFTER;
192     double obj[DROP];
193     double kbad = 0, kgood = 0;
194     if (strategy & 128) after = 999999; /* no acceleration at all */
195     for (i = 0; i < DROP; i++) {
196          obj[i] = 1.0e70;
197     }
198     //#define TWO_GOES
199#ifdef TWO_GOES
200     double * pi2 = new double [nrows];
201     double * rowsol2 = new double [nrows];
202#endif
203     allsum = new double * [nsolve];
204     aX = new double * [nsolve];
205     aworkX = new double * [nsolve];
206     thetaX = new double[nsolve];
207     vX = new double[nsolve];
208     bX = new double[nsolve];
209     djX = new double[nsolve];
210     allsum[0] = pi;
211     for (i = 0; i < nsolve; i++) {
212          if (i) allsum[i] = new double[nrows];
213          aX[i] = new double[nsolve];
214          aworkX[i] = new double[nsolve];
215     }
216     /* check = rows */
217     for (i = 0; i < nrows; i++) {
218          if (rowupper[i] - rowlower[i] > tolerance) {
219               extraBlock++;
220          }
221     }
222     cost = new double[ncols];
223     memset(rowsol, 0, nrows * sizeof(double));
224     for (i = 0; i < ncols; i++) {
225          CoinBigIndex j;
226          double value = origcost[i] * maxmin;
227          double value2 = colsol[i];
228          if (elemnt) {
229               for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
230                    int irow = row[j];
231                    value += elemnt[j] * lambda[irow];
232                    rowsol[irow] += elemnt[j] * value2;
233               }
234          } else {
235               for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
236                    int irow = row[j];
237                    value += lambda[irow];
238                    rowsol[irow] += value2;
239               }
240          }
241          cost[i] = value;
242     }
243     if (extraBlock) {
244          rowExtra = new int[extraBlock];
245          solExtra = new double[extraBlock];
246          elemExtra = new double[extraBlock];
247          upperExtra = new double[extraBlock];
248          costExtra = new double[extraBlock];
249          saveExtra = new double[extraBlock];
250          extraBlock = 0;
251          int nbad = 0;
252          for (i = 0; i < nrows; i++) {
253               if (rowupper[i] - rowlower[i] > tolerance) {
254                    double smaller, difference;
255                    double value;
256                    saveExtra[extraBlock] = rowupper[i];
257                    if (fabs(rowupper[i]) > fabs(rowlower[i])) {
258                         smaller = rowlower[i];
259                         value = -1.0;
260                    } else {
261                         smaller = rowupper[i];
262                         value = 1.0;
263                    }
264                    if (fabs(smaller) > 1.0e10) {
265                         if (!nbad)
266                              printf("Can't handle rows where both bounds >1.0e10 %d %g\n",
267                                     i, smaller);
268                         nbad++;
269                         if (rowupper[i] < 0.0 || rowlower[i] > 0.0)
270                              abort();
271                         if (fabs(rowupper[i]) > fabs(rowlower[i])) {
272                              rowlower[i] = -0.9e10;
273                              smaller = rowlower[i];
274                              value = -1.0;
275                         } else {
276                              rowupper[i] = 0.9e10;
277                              saveExtra[extraBlock] = rowupper[i];
278                              smaller = rowupper[i];
279                              value = 1.0;
280                         }
281                    }
282                    difference = rowupper[i] - rowlower[i];
283                    difference = CoinMin(difference, 1.0e31);
284                    rowupper[i] = smaller;
285                    elemExtra[extraBlock] = value;
286                    solExtra[extraBlock] = (rowupper[i] - rowsol[i]) / value;
287                    if (solExtra[extraBlock] < 0.0) solExtra[extraBlock] = 0.0;
288                    if (solExtra[extraBlock] > difference) solExtra[extraBlock] = difference;
289                    costExtra[extraBlock] = lambda[i] * value;
290                    upperExtra[extraBlock] = difference;
291                    rowsol[i] += value * solExtra[extraBlock];
292                    rowExtra[extraBlock++] = i;
293               }
294          }
295          if (nbad)
296               printf("%d bad values - results may be wrong\n", nbad);
297     }
298     for (i = 0; i < nrows; i++) {
299          offset += lambda[i] * rowsol[i];
300     }
301     if ((strategy & 256) != 0) {
302          /* save best solution */
303          saveSol = new double[ncols];
304          CoinMemcpyN(colsol, ncols, saveSol);
305          if (extraBlock) {
306               useCostExtra = new double[extraBlock];
307               memset(useCostExtra, 0, extraBlock * sizeof(double));
308          }
309          useCost = origcost;
310          useOffset = saveOffset;
311     } else {
312          useCostExtra = costExtra;
313          useCost = cost;
314          useOffset = offset;
315     }
316     ncolx = ncols + extraBlock;
317     for (i = 0; i < HISTORY + 1; i++) {
318          history[i] = new double[ncolx];
319     }
320     for (i = 0; i < DJTEST; i++) {
321          djSave[i] = 1.0e30;
322     }
323     for (i = 0; i < ncols; i++) {
324          if (upper[i] - lower[i]) {
325               statusSave[i] = 0;
326          } else {
327               statusSave[i] = 1;
328          }
329     }
330     // for two pass method
331     int start[2];
332     int stop[2];
333     int direction = -1;
334     start[0] = 0;
335     stop[0] = ncols;
336     start[1] = 0;
337     stop[1] = 0;
338     iter = 0;
339     for (; iter < maxIts; iter++) {
340          double sum1 = 0.0, sum2 = 0.0, djtot;
341          double lastObj = 1.0e70;
342          int good = 0, doScale = 0;
343          if (strategy & 16) {
344               int ii = iter / EVERY + 1;
345               ii = ii * EVERY;
346               if (iter > ii - HISTORY * 2 && (iter & 1) == 0) {
347                    double * x = history[HISTORY-1];
348                    for (i = HISTORY - 1; i > 0; i--) {
349                         history[i] = history[i-1];
350                    }
351                    history[0] = x;
352                    CoinMemcpyN(colsol, ncols, history[0]);
353                    CoinMemcpyN(solExtra, extraBlock, history[0] + ncols);
354               }
355          }
356          if ((iter % SAVEHISTORY) == 0 || doFull) {
357               if ((strategy & 16) == 0) {
358                    double * x = history[HISTORY-1];
359                    for (i = HISTORY - 1; i > 0; i--) {
360                         history[i] = history[i-1];
361                    }
362                    history[0] = x;
363                    CoinMemcpyN(colsol, ncols, history[0]);
364                    CoinMemcpyN(solExtra, extraBlock, history[0] + ncols);
365               }
366          }
367          /* start full try */
368          if ((iter % EVERY) == 0 || doFull) {
369               // for next pass
370               direction = - direction;
371               // randomize.
372               // The cast is to avoid gcc compiler warning
373               int kcol = static_cast<int>(ncols * randomNumberGenerator->randomDouble());
374               if (kcol == ncols)
375                    kcol = ncols - 1;
376               if (direction > 0) {
377                    start[0] = kcol;
378                    stop[0] = ncols;
379                    start[1] = 0;
380                    stop[1] = kcol;
381               } else {
382                    start[0] = kcol;
383                    stop[0] = -1;
384                    start[1] = ncols - 1;
385                    stop[1] = kcol;
386               }
387               int itry = 0;
388               /*if ((strategy&16)==0) {
389                double * x=history[HISTORY-1];
390                for (i=HISTORY-1;i>0;i--) {
391                history[i]=history[i-1];
392                }
393                history[0]=x;
394                CoinMemcpyN(colsol,ncols,history[0]);
395                 CoinMemcpyN(solExtra,extraBlock,history[0]+ncols);
396                }*/
397               while (!good) {
398                    itry++;
399#define MAXTRY 5
400                    if (iter > after && doScale < 2 && itry < MAXTRY) {
401                         /* now full one */
402                         for (i = 0; i < nrows; i++) {
403                              rowsol[i] = -rowupper[i];
404                         }
405                         sum2 = 0.0;
406                         objvalue = 0.0;
407                         memset(pi, 0, nrows * sizeof(double));
408                         djtot = 0.0;
409                         {
410                              double * theta = thetaX;
411                              double * dj = djX;
412                              double * b = bX;
413                              double ** a = aX;
414                              double ** awork = aworkX;
415                              double * v = vX;
416                              double c;
417#ifdef FIT
418                              int ntot = 0, nsign = 0, ngood = 0, mgood[4] = {0, 0, 0, 0};
419                              double diff1, diff2, val0, val1, val2, newValue;
420                              CoinMemcpyN(colsol, ncols, history[HISTORY-1]);
421                              CoinMemcpyN(solExtra, extraBlock, history[HISTORY-1] + ncols);
422#endif
423                              dj[0] = 0.0;
424                              for (i = 1; i < nsolve; i++) {
425                                   dj[i] = 0.0;
426                                   memset(allsum[i], 0, nrows * sizeof(double));
427                              }
428                              for (i = 0; i < ncols; i++) {
429                                   double value2 = colsol[i];
430                                   if (value2 > lower[i] + tolerance) {
431                                        if(value2 < (upper[i] - tolerance)) {
432                                             int k;
433                                             objvalue += value2 * cost[i];
434#ifdef FIT
435                                             ntot++;
436                                             val0 = history[0][i];
437                                             val1 = history[1][i];
438                                             val2 = history[2][i];
439                                             diff1 = val0 - val1;
440                                             diff2 = val1 - val2;
441                                             if (diff1*diff2 >= 0.0) {
442                                                  nsign++;
443                                                  if (fabs(diff1) < fabs(diff2)) {
444                                                       int ii = static_cast<int>(fabs(4.0 * diff1 / diff2));
445                                                       if (ii == 4) ii = 3;
446                                                       mgood[ii]++;
447                                                       ngood++;
448                                                  }
449                                                  if (fabs(diff1) < 0.75 * fabs(diff2)) {
450                                                       newValue = val1 + (diff1 * diff2) / (diff2 - diff1);
451                                                  } else {
452                                                       newValue = val1 + 4.0 * diff1;
453                                                  }
454                                             } else {
455                                                  newValue = 0.333333333 * (val0 + val1 + val2);
456                                             }
457                                             if (newValue > upper[i] - tolerance) {
458                                                  newValue = upper[i];
459                                             } else if (newValue < lower[i] + tolerance) {
460                                                  newValue = lower[i];
461                                             }
462                                             history[HISTORY-1][i] = newValue;
463#endif
464                                             for (k = 0; k < HISTORY - 1; k++) {
465                                                  value = history[k][i] - history[k+1][i];
466                                                  dj[k] += value * cost[i];
467                                                  v[k] = value;
468                                             }
469                                             if (elemnt) {
470                                                  for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
471                                                       int irow = row[j];
472                                                       for (k = 0; k < HISTORY - 1; k++) {
473                                                            allsum[k][irow] += elemnt[j] * v[k];
474                                                       }
475                                                       rowsol[irow] += elemnt[j] * value2;
476                                                  }
477                                             } else {
478                                                  for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
479                                                       int irow = row[j];
480                                                       for (k = 0; k < HISTORY - 1; k++) {
481                                                            allsum[k][irow] += v[k];
482                                                       }
483                                                       rowsol[irow] += value2;
484                                                  }
485                                             }
486                                        } else {
487                                             /* at ub */
488                                             colsol[i] = upper[i];
489                                             value2 = colsol[i];
490                                             objvalue += value2 * cost[i];
491                                             if (elemnt) {
492                                                  for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
493                                                       int irow = row[j];
494                                                       rowsol[irow] += elemnt[j] * value2;
495                                                  }
496                                             } else {
497                                                  for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
498                                                       int irow = row[j];
499                                                       rowsol[irow] += value2;
500                                                  }
501                                             }
502                                        }
503                                   } else {
504                                        /* at lb */
505                                        if (value2) {
506                                             objvalue += value2 * cost[i];
507                                             if (elemnt) {
508                                                  for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
509                                                       int irow = row[j];
510                                                       rowsol[irow] += elemnt[j] * value2;
511                                                  }
512                                             } else {
513                                                  for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
514                                                       int irow = row[j];
515                                                       rowsol[irow] += value2;
516                                                  }
517                                             }
518                                        }
519                                   }
520                              }
521#ifdef FIT
522                              /*printf("total %d, same sign %d, good %d %d %d %d %d\n",
523                                ntot,nsign,ngood,mgood[0],mgood[1],mgood[2],mgood[3]);*/
524#endif
525                              if (extraBlock) {
526                                   for (i = 0; i < extraBlock; i++) {
527                                        double element = elemExtra[i];
528                                        int irow = rowExtra[i];
529                                        objvalue += solExtra[i] * costExtra[i];
530                                        if (solExtra[i] > tolerance
531                                                  && solExtra[i] < (upperExtra[i] - tolerance)) {
532                                             double value2 = solExtra[i];
533                                             int k;
534                                             for (k = 0; k < HISTORY - 1; k++) {
535                                                  value = history[k][i+ncols] - history[k+1][i+ncols];
536                                                  dj[k] += value * costExtra[i];
537                                                  allsum[k][irow] += element * value;
538                                             }
539                                             rowsol[irow] += element * value2;
540                                        } else {
541                                             double value2 = solExtra[i];
542                                             double element = elemExtra[i];
543                                             int irow = rowExtra[i];
544                                             rowsol[irow] += element * value2;
545                                        }
546                                   }
547                              }
548#ifdef NULLVECTOR
549                              if ((strategy & 64)) {
550                                   double djVal = dj[0];
551                                   for (i = 0; i < ncols - nrows; i++) {
552                                        double value2 = colsol[i];
553                                        if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance) {
554                                             value = history[0][i] - history[1][i];
555                                        } else {
556                                             value = 0.0;
557                                        }
558                                        history[HISTORY][i] = value;
559                                   }
560                                   for (; i < ncols; i++) {
561                                        double value2 = colsol[i];
562                                        double delta;
563                                        int irow = i - (ncols - nrows);
564                                        double oldSum = allsum[0][irow];;
565                                        if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance) {
566                                             delta = history[0][i] - history[1][i];
567                                        } else {
568                                             delta = 0.0;
569                                        }
570                                        djVal -= delta * cost[i];
571                                        oldSum -= delta;
572                                        delta = - oldSum;
573                                        djVal += delta * cost[i];
574                                        history[HISTORY][i] = delta;
575                                   }
576                                   dj[HISTORY-1] = djVal;
577                                   djVal = 0.0;
578                                   for (i = 0; i < ncols; i++) {
579                                        double value2 = colsol[i];
580                                        if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance ||
581                                                  i >= ncols - nrows) {
582                                             int k;
583                                             value = history[HISTORY][i];
584                                             djVal += value * cost[i];
585                                             for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
586                                                  int irow = row[j];
587                                                  allsum[nsolve-1][irow] += value;
588                                             }
589                                        }
590                                   }
591                                   printf("djs %g %g\n", dj[HISTORY-1], djVal);
592                              }
593#endif
594                              for (i = 0; i < nsolve; i++) {
595                                   int j;
596                                   b[i] = 0.0;
597                                   for (j = 0; j < nsolve; j++) {
598                                        a[i][j] = 0.0;
599                                   }
600                              }
601                              c = 0.0;
602                              for (i = 0; i < nrows; i++) {
603                                   double value = rowsol[i];
604                                   for (k = 0; k < nsolve; k++) {
605                                        v[k] = allsum[k][i];
606                                        b[k] += v[k] * value;
607                                   }
608                                   c += value * value;
609                                   for (k = 0; k < nsolve; k++) {
610                                        for (j = k; j < nsolve; j++) {
611                                             a[k][j] += v[k] * v[j];
612                                        }
613                                   }
614                              }
615                              sum2 = c;
616                              if (itry == 1) {
617                                   lastObj = objvalue + weight * sum2;
618                              }
619                              for (k = 0; k < nsolve; k++) {
620                                   b[k] = - (weight * b[k] + 0.5 * dj[k]);
621                                   for (j = k; j < nsolve; j++) {
622                                        a[k][j] *= weight;
623                                        a[j][k] = a[k][j];
624                                   }
625                              }
626                              c *= weight;
627                              for (k = 0; k < nsolve; k++) {
628                                   theta[k] = b[k];
629                              }
630                              solveSmall(nsolve, a, awork, theta);
631                              if ((strategy & 64) != 0) {
632                                   value = 10.0;
633                                   for (k = 0; k < nsolve; k++) {
634                                        value  = CoinMax(value, fabs(theta[k]));
635                                   }
636                                   if (value > 10.0 && ((logLevel_ & 4) != 0)) {
637                                        printf("theta %g %g %g\n", theta[0], theta[1], theta[2]);
638                                   }
639                                   value = 10.0 / value;
640                                   for (k = 0; k < nsolve; k++) {
641                                        theta[k] *= value;
642                                   }
643                              }
644                              for (i = 0; i < ncolx; i++) {
645                                   double valueh = 0.0;
646                                   for (k = 0; k < HISTORY - 1; k++) {
647                                        value = history[k][i] - history[k+1][i];
648                                        valueh += value * theta[k];
649                                   }
650#ifdef NULLVECTOR
651                                   value = history[HISTORY][i];
652                                   valueh += value * theta[HISTORY-1];
653#endif
654                                   history[HISTORY][i] = valueh;
655                              }
656                         }
657#ifdef NULLVECTOR
658                         if ((strategy & 64)) {
659                              for (i = 0; i < ncols - nrows; i++) {
660                                   if (colsol[i] <= lower[i] + tolerance
661                                             || colsol[i] >= (upper[i] - tolerance)) {
662                                        history[HISTORY][i] = 0.0;;
663                                   }
664                              }
665                              tolerance = -tolerance; /* switch off test */
666                         }
667#endif
668                         if (!doScale) {
669                              for (i = 0; i < ncols; i++) {
670                                   if (colsol[i] > lower[i] + tolerance
671                                             && colsol[i] < (upper[i] - tolerance)) {
672                                        value = history[HISTORY][i];
673                                        colsol[i] += value;
674                                        if (colsol[i] < lower[i] + tolerance) {
675                                             colsol[i] = lower[i];
676                                        } else if (colsol[i] > upper[i] - tolerance) {
677                                             colsol[i] = upper[i];
678                                        }
679                                   }
680                              }
681                              if (extraBlock) {
682                                   for (i = 0; i < extraBlock; i++) {
683                                        if (solExtra[i] > tolerance
684                                                  && solExtra[i] < (upperExtra[i] - tolerance)) {
685                                             value = history[HISTORY][i+ncols];
686                                             solExtra[i] += value;
687                                             if (solExtra[i] < 0.0) {
688                                                  solExtra[i] = 0.0;
689                                             } else if (solExtra[i] > upperExtra[i]) {
690                                                  solExtra[i] = upperExtra[i];
691                                             }
692                                        }
693                                   }
694                              }
695                         } else {
696                              double theta = 1.0;
697                              double saveTheta = theta;
698                              for (i = 0; i < ncols; i++) {
699                                   if (colsol[i] > lower[i] + tolerance
700                                             && colsol[i] < (upper[i] - tolerance)) {
701                                        value = history[HISTORY][i];
702                                        if (value > 0) {
703                                             if (theta * value + colsol[i] > upper[i]) {
704                                                  theta = (upper[i] - colsol[i]) / value;
705                                             }
706                                        } else if (value < 0) {
707                                             if (colsol[i] + theta * value < lower[i]) {
708                                                  theta = (lower[i] - colsol[i]) / value;
709                                             }
710                                        }
711                                   }
712                              }
713                              if (extraBlock) {
714                                   for (i = 0; i < extraBlock; i++) {
715                                        if (solExtra[i] > tolerance
716                                                  && solExtra[i] < (upperExtra[i] - tolerance)) {
717                                             value = history[HISTORY][i+ncols];
718                                             if (value > 0) {
719                                                  if (theta * value + solExtra[i] > upperExtra[i]) {
720                                                       theta = (upperExtra[i] - solExtra[i]) / value;
721                                                  }
722                                             } else if (value < 0) {
723                                                  if (solExtra[i] + theta * value < 0.0) {
724                                                       theta = -solExtra[i] / value;
725                                                  }
726                                             }
727                                        }
728                                   }
729                              }
730                              if ((iter % 100 == 0) && (logLevel_ & 8) != 0) {
731                                   if (theta < saveTheta) {
732                                        printf(" - modified theta %g\n", theta);
733                                   }
734                              }
735                              for (i = 0; i < ncols; i++) {
736                                   if (colsol[i] > lower[i] + tolerance
737                                             && colsol[i] < (upper[i] - tolerance)) {
738                                        value = history[HISTORY][i];
739                                        colsol[i] += value * theta;
740                                   }
741                              }
742                              if (extraBlock) {
743                                   for (i = 0; i < extraBlock; i++) {
744                                        if (solExtra[i] > tolerance
745                                                  && solExtra[i] < (upperExtra[i] - tolerance)) {
746                                             value = history[HISTORY][i+ncols];
747                                             solExtra[i] += value * theta;
748                                        }
749                                   }
750                              }
751                         }
752#ifdef NULLVECTOR
753                         tolerance = fabs(tolerance); /* switch back on */
754#endif
755                         if ((iter % 100) == 0 && (logLevel_ & 8) != 0) {
756                              printf("\n");
757                         }
758                    }
759                    good = 1;
760                    result = objval(nrows, ncols, rowsol, colsol, pi, djs, useCost,
761                                    rowlower, rowupper, lower, upper,
762                                    elemnt, row, columnStart, length, extraBlock, rowExtra,
763                                    solExtra, elemExtra, upperExtra, useCostExtra,
764                                    weight);
765                    weightedObj = result.weighted;
766                    if (!iter) saveValue = weightedObj;
767                    objvalue = result.objval;
768                    sum1 = result.infeas;
769                    if (saveSol) {
770                         if (result.weighted < bestSol) {
771                              printf("%d %g better than %g\n", iter,
772                                     result.weighted * maxmin - useOffset, bestSol * maxmin - useOffset);
773                              bestSol = result.weighted;
774                              CoinMemcpyN(colsol, ncols, saveSol);
775                         }
776                    }
777#ifdef FITz
778                    if (iter > after) {
779                         IdiotResult result2;
780                         double ww, oo, ss;
781                         if (extraBlock) abort();
782                         result2 = objval(nrows, ncols, row2, sol2, pi2, djs, cost,
783                                          rowlower, rowupper, lower, upper,
784                                          elemnt, row, columnStart, extraBlock, rowExtra,
785                                          solExtra, elemExtra, upperExtra, costExtra,
786                                          weight);
787                         ww = result2.weighted;
788                         oo = result2.objval;
789                         ss = result2.infeas;
790                         printf("wobj %g obj %g inf %g last %g\n", ww, oo, ss, lastObj);
791                         if (ww < weightedObj && ww < lastObj) {
792                              printf(" taken");
793                              ntaken++;
794                              saving += weightedObj - ww;
795                              weightedObj = ww;
796                              objvalue = oo;
797                              sum1 = ss;
798                              CoinMemcpyN(row2, nrows, rowsol);
799                              CoinMemcpyN(pi2, nrows, pi);
800                              CoinMemcpyN(sol2, ncols, colsol);
801                              result = objval(nrows, ncols, rowsol, colsol, pi, djs, cost,
802                                              rowlower, rowupper, lower, upper,
803                                              elemnt, row, columnStart, extraBlock, rowExtra,
804                                              solExtra, elemExtra, upperExtra, costExtra,
805                                              weight);
806                              weightedObj = result.weighted;
807                              objvalue = result.objval;
808                              sum1 = result.infeas;
809                              if (ww < weightedObj) abort();
810                         } else {
811                              printf(" not taken");
812                              nottaken++;
813                         }
814                    }
815#endif
816                    /*printf("%d %g %g %g %g\n",itry,lastObj,weightedObj,objvalue,sum1);*/
817                    if (weightedObj > lastObj + 1.0e-4 && itry < MAXTRY) {
818                         if((logLevel_ & 16) != 0 && doScale) {
819                              printf("Weighted objective from %g to %g **** bad move\n",
820                                     lastObj, weightedObj);
821                         }
822                         if (doScale) {
823                              good = 1;
824                         }
825                         if ((strategy & 3) == 1) {
826                              good = 0;
827                              if (weightedObj > lastObj + djExit) {
828                                   if ((logLevel_ & 16) != 0) {
829                                        printf("Weighted objective from %g to %g ?\n", lastObj, weightedObj);
830                                   }
831                                   CoinMemcpyN(history[0], ncols, colsol);
832                                   CoinMemcpyN(history[0] + ncols, extraBlock, solExtra);
833                                   good = 1;
834                              }
835                         } else if ((strategy & 3) == 2) {
836                              if (weightedObj > lastObj + 0.1 * maxDj) {
837                                   CoinMemcpyN(history[0], ncols, colsol);
838                                   CoinMemcpyN(history[0] + ncols, extraBlock, solExtra);
839                                   doScale++;
840                                   good = 0;
841                              }
842                         } else if ((strategy & 3) == 3) {
843                              if (weightedObj > lastObj + 0.001 * maxDj) {
844                                   /*doScale++;*/
845                                   good = 0;
846                              }
847                         }
848                    }
849               }
850               if ((iter % checkFrequency_) == 0) {
851                    double best = weightedObj;
852                    double test = obj[0];
853                    for (i = 1; i < DROP; i++) {
854                         obj[i-1] = obj[i];
855                         if (best > obj[i]) best = obj[i];
856                    }
857                    obj[DROP-1] = best;
858                    if (test - best < drop && (strategy & 8) == 0) {
859                         if ((logLevel_ & 8) != 0) {
860                              printf("Exiting as drop in %d its is %g after %d iterations\n",
861                                     DROP * checkFrequency_, test - best, iter);
862                         }
863                         goto RETURN;
864                    }
865               }
866               if ((iter % logFreq_) == 0) {
867                    double piSum = 0.0;
868                    for (i = 0; i < nrows; i++) {
869                         piSum += (rowsol[i] + rowupper[i]) * pi[i];
870                    }
871                    if ((logLevel_ & 2) != 0) {
872                         printf("%d Infeas %g, obj %g - wtObj %g dual %g maxDj %g\n",
873                                iter, sum1, objvalue * maxmin - useOffset, weightedObj - useOffset,
874                                piSum * maxmin - useOffset, maxDj);
875                    }
876               }
877               CoinMemcpyN(statusSave, ncols, statusWork);
878               nflagged = 0;
879          }
880          nChange = 0;
881          doFull = 0;
882          maxDj = 0.0;
883          // go through forwards or backwards and starting at odd places
884          int itry;
885#ifdef TWO_GOES
886          memcpy(pi2, pi, nrows * sizeof(double));
887          memcpy(rowsol2, rowsol, nrows * sizeof(double));
888#endif
889          for (itry = 0; itry < 2; itry++) {
890               int icol = start[itry];
891               int istop = stop[itry];
892#ifdef TWO_GOES
893               for (int iPar = 0; iPar < 2; iPar++) {
894                    double * temp = pi;
895                    pi = pi2;
896                    pi2 = temp;
897                    temp = rowsol;
898                    rowsol = rowsol2;
899                    rowsol2 = temp;
900                    if (iPar == 0) {
901                         istop = (icol + istop) >> 1;
902                    } else {
903                         icol = istop;
904                         istop = stop[itry];
905                    }
906#endif
907                    for (; icol != istop; icol += direction) {
908                         if (!statusWork[icol]) {
909                              CoinBigIndex j;
910                              double value = colsol[icol];
911                              double djval = cost[icol];
912                              double djval2, value2;
913                              double theta, a, b, c;
914                              if (elemnt) {
915                                   for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
916                                        int irow = row[j];
917                                        djval -= elemnt[j] * pi[irow];
918                                   }
919                              } else {
920                                   for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
921                                        int irow = row[j];
922                                        djval -= pi[irow];
923                                   }
924                              }
925                              /*printf("xx iter %d seq %d djval %g value %g\n",
926                                iter,i,djval,value);*/
927                              if (djval > 1.0e-5) {
928                                   value2 = (lower[icol] - value);
929                              } else {
930                                   value2 = (upper[icol] - value);
931                              }
932                              djval2 = djval * value2;
933                              djval = fabs(djval);
934                              if (djval > djTol) {
935                                   if (djval2 < -1.0e-4) {
936                                        double valuep, thetap;
937                                        nChange++;
938                                        if (djval > maxDj) maxDj = djval;
939                                        /*if (djval>3.55e6) {
940                                                        printf("big\n");
941                                                        }*/
942                                        a = 0.0;
943                                        b = 0.0;
944                                        c = 0.0;
945                                        djval2 = cost[icol];
946                                        if (elemnt) {
947                                             for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
948                                                  int irow = row[j];
949                                                  double value = rowsol[irow];
950                                                  c += value * value;
951                                                  a += elemnt[j] * elemnt[j];
952                                                  b += value * elemnt[j];
953                                             }
954                                        } else {
955                                             for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
956                                                  int irow = row[j];
957                                                  double value = rowsol[irow];
958                                                  c += value * value;
959                                                  a += 1.0;
960                                                  b += value;
961                                             }
962                                        }
963                                        a *= weight;
964                                        b = b * weight + 0.5 * djval2;
965                                        c *= weight;
966                                        /* solve */
967                                        theta = -b / a;
968                                        if ((strategy & 4) != 0) {
969                                             value2 = a * theta * theta + 2.0 * b * theta;
970                                             thetap = 2.0 * theta;
971                                             valuep = a * thetap * thetap + 2.0 * b * thetap;
972                                             if (valuep < value2 + djTol) {
973                                                  theta = thetap;
974                                                  kgood++;
975                                             } else {
976                                                  kbad++;
977                                             }
978                                        }
979                                        if (theta > 0.0) {
980                                             if (theta < upper[icol] - colsol[icol]) {
981                                                  value2 = theta;
982                                             } else {
983                                                  value2 = upper[icol] - colsol[icol];
984                                             }
985                                        } else {
986                                             if (theta > lower[icol] - colsol[icol]) {
987                                                  value2 = theta;
988                                             } else {
989                                                  value2 = lower[icol] - colsol[icol];
990                                             }
991                                        }
992                                        colsol[icol] += value2;
993                                        objvalue += cost[icol] * value2;
994                                        if (elemnt) {
995                                             for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
996                                                  int irow = row[j];
997                                                  double value;
998                                                  rowsol[irow] += elemnt[j] * value2;
999                                                  value = rowsol[irow];
1000                                                  pi[irow] = -2.0 * weight * value;
1001                                             }
1002                                        } else {
1003                                             for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
1004                                                  int irow = row[j];
1005                                                  double value;
1006                                                  rowsol[irow] += value2;
1007                                                  value = rowsol[irow];
1008                                                  pi[irow] = -2.0 * weight * value;
1009                                             }
1010                                        }
1011                                   } else {
1012                                        /* dj but at bound */
1013                                        if (djval > djFlag) {
1014                                             statusWork[icol] = 1;
1015                                             nflagged++;
1016                                        }
1017                                   }
1018                              }
1019                         }
1020                    }
1021#ifdef TWO_GOES
1022               }
1023#endif
1024          }
1025#ifdef TWO_GOES
1026          for (i = 0; i < nrows; i++) {
1027               rowsol[i] = 0.5 * (rowsol[i] + rowsol2[i]);
1028               pi[i] = 0.5 * (pi[i] + pi2[i]);
1029          }
1030#endif
1031          if (extraBlock) {
1032               for (i = 0; i < extraBlock; i++) {
1033                    double value = solExtra[i];
1034                    double djval = costExtra[i];
1035                    double djval2, value2;
1036                    double element = elemExtra[i];
1037                    double theta, a, b, c;
1038                    int irow = rowExtra[i];
1039                    djval -= element * pi[irow];
1040                    /*printf("xxx iter %d extra %d djval %g value %g\n",
1041                          iter,irow,djval,value);*/
1042                    if (djval > 1.0e-5) {
1043                         value2 = -value;
1044                    } else {
1045                         value2 = (upperExtra[i] - value);
1046                    }
1047                    djval2 = djval * value2;
1048                    if (djval2 < -1.0e-4 && fabs(djval) > djTol) {
1049                         nChange++;
1050                         a = 0.0;
1051                         b = 0.0;
1052                         c = 0.0;
1053                         djval2 = costExtra[i];
1054                         value = rowsol[irow];
1055                         c += value * value;
1056                         a += element * element;
1057                         b += element * value;
1058                         a *= weight;
1059                         b = b * weight + 0.5 * djval2;
1060                         c *= weight;
1061                         /* solve */
1062                         theta = -b / a;
1063                         if (theta > 0.0) {
1064                              value2 = CoinMin(theta, upperExtra[i] - solExtra[i]);
1065                         } else {
1066                              value2 = CoinMax(theta, -solExtra[i]);
1067                         }
1068                         solExtra[i] += value2;
1069                         rowsol[irow] += element * value2;
1070                         value = rowsol[irow];
1071                         pi[irow] = -2.0 * weight * value;
1072                    }
1073               }
1074          }
1075          if ((iter % 10) == 2) {
1076               for (i = DJTEST - 1; i > 0; i--) {
1077                    djSave[i] = djSave[i-1];
1078               }
1079               djSave[0] = maxDj;
1080               largestDj = CoinMax(largestDj, maxDj);
1081               smallestDj = CoinMin(smallestDj, maxDj);
1082               for (i = DJTEST - 1; i > 0; i--) {
1083                    maxDj += djSave[i];
1084               }
1085               maxDj = maxDj / static_cast<double> (DJTEST);
1086               if (maxDj < djExit && iter > 50) {
1087                    //printf("Exiting on low dj %g after %d iterations\n",maxDj,iter);
1088                    break;
1089               }
1090               if (nChange < 100) {
1091                    djTol *= 0.5;
1092               }
1093          }
1094     }
1095RETURN:
1096     if (kgood || kbad) {
1097          printf("%g good %g bad\n", kgood, kbad);
1098     }
1099     result = objval(nrows, ncols, rowsol, colsol, pi, djs, useCost,
1100                     rowlower, rowupper, lower, upper,
1101                     elemnt, row, columnStart, length, extraBlock, rowExtra,
1102                     solExtra, elemExtra, upperExtra, useCostExtra,
1103                     weight);
1104     result.djAtBeginning = largestDj;
1105     result.djAtEnd = smallestDj;
1106     result.dropThis = saveValue - result.weighted;
1107     if (saveSol) {
1108          if (result.weighted < bestSol) {
1109               bestSol = result.weighted;
1110               CoinMemcpyN(colsol, ncols, saveSol);
1111          } else {
1112               printf("restoring previous - now %g best %g\n",
1113                      result.weighted * maxmin - useOffset, bestSol * maxmin - useOffset);
1114          }
1115     }
1116     if (saveSol) {
1117          if (extraBlock) {
1118               delete [] useCostExtra;
1119          }
1120          CoinMemcpyN(saveSol, ncols, colsol);
1121          delete [] saveSol;
1122     }
1123     for (i = 0; i < nsolve; i++) {
1124          if (i) delete [] allsum[i];
1125          delete [] aX[i];
1126          delete [] aworkX[i];
1127     }
1128     delete [] thetaX;
1129     delete [] djX;
1130     delete [] bX;
1131     delete [] vX;
1132     delete [] aX;
1133     delete [] aworkX;
1134     delete [] allsum;
1135     delete [] cost;
1136#ifdef TWO_GOES
1137     delete [] pi2 ;
1138     delete [] rowsol2 ;
1139#endif
1140     for (i = 0; i < HISTORY + 1; i++) {
1141          delete [] history[i];
1142     }
1143     delete [] statusSave;
1144     /* do original costs objvalue*/
1145     result.objval = 0.0;
1146     for (i = 0; i < ncols; i++) {
1147          result.objval += colsol[i] * origcost[i];
1148     }
1149     if (extraBlock) {
1150          for (i = 0; i < extraBlock; i++) {
1151               int irow = rowExtra[i];
1152               rowupper[irow] = saveExtra[i];
1153          }
1154          delete [] rowExtra;
1155          delete [] solExtra;
1156          delete [] elemExtra;
1157          delete [] upperExtra;
1158          delete [] costExtra;
1159          delete [] saveExtra;
1160     }
1161     result.iteration = iter;
1162     result.objval -= saveOffset;
1163     result.weighted = result.objval + weight * result.sumSquared;
1164     return result;
1165}
Note: See TracBrowser for help on using the repository browser.