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

Last change on this file since 754 was 754, checked in by andreasw, 14 years ago

first version

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