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

Last change on this file since 2470 was 2385, checked in by unxusr, 10 months ago

formatting

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