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

Last change on this file since 2030 was 2030, checked in by forrest, 6 years ago

fix some ampl stuff, add ClpSolver? and a few fixes

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