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

Last change on this file since 1141 was 1141, checked in by forrest, 12 years ago

for threaded random numbers

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