Changeset 102 for trunk


Ignore:
Timestamp:
Jan 22, 2003 3:42:41 PM (17 years ago)
Author:
forrest
Message:

Crash

Location:
trunk
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpMessage.cpp

    r93 r102  
    7171  {CLP_IMPORT_ERRORS,3001,1," There were %d errors when importing model from %s"},
    7272  {CLP_EMPTY_PROBLEM,3002,0,"Not solving empty problem - %d rows, %d columns and %d elements"},
     73  {CLP_CRASH,28,1,"Crash put %d variables in basis, %d dual infeasibilities"},
    7374  {CLP_DUMMY_END,999999,0,""}
    7475};
  • trunk/ClpSimplex.cpp

    r98 r102  
    32423242ClpSimplex::crash(double gap,int pivot)
    32433243{
    3244   assert(!pivot); // rest not coded yet
    32453244  assert(!rowObjective_); // not coded
    32463245  int iColumn;
    32473246  int numberBad=0;
    32483247  int numberBasic=0;
     3248  double dualTolerance=dblParam_[ClpDualTolerance];
     3249  //double primalTolerance=dblParam_[ClpPrimalTolerance];
    32493250
    32503251  // If no basis then make all slack one
     
    32613262  } else {
    32623263    double * dj = new double [numberColumns_];
    3263     double objectiveValue=0.0;
     3264    double * solution = columnActivity_;
     3265    //double objectiveValue=0.0;
    32643266    int iColumn;
    32653267    for (iColumn=0;iColumn<numberColumns_;iColumn++)
     
    32743276          atLower=false;
    32753277          setColumnStatus(iColumn,atUpperBound);
     3278          solution[iColumn]=upperBound;
    32763279        } else {
    32773280          atLower=true;
    32783281          setColumnStatus(iColumn,atLowerBound);
     3282          solution[iColumn]=lowerBound;
    32793283        }
    32803284        if (dj[iColumn]<0.0) {
     
    32853289              columnActivity_[iColumn]=upperBound;
    32863290              setColumnStatus(iColumn,atUpperBound);
    3287             } else if (dj[iColumn]<-dualTolerance_) {
     3291            } else if (dj[iColumn]<-dualTolerance) {
    32883292              numberBad++;
    32893293            }
     
    32963300              columnActivity_[iColumn]=lowerBound;
    32973301              setColumnStatus(iColumn,atLowerBound);
    3298             } else if (dj[iColumn]>dualTolerance_) {
     3302            } else if (dj[iColumn]>dualTolerance) {
    32993303              numberBad++;
    33003304            }
     
    33043308        // free
    33053309        setColumnStatus(iColumn,isFree);
    3306         if (fabs(dj[iColumn])>dualTolerance_)
     3310        if (fabs(dj[iColumn])>dualTolerance)
    33073311          numberBad++;
    33083312      }
    33093313    }
    3310     if (numberBad) {
     3314    if (numberBad||pivot) {
    33113315      if (!pivot) {
    33123316        delete [] dj;
    33133317        return 1;
    33143318      } else {
    3315 #if 0
    33163319        // see if can be made dual feasible with gubs etc
    33173320        double * pi = new double[numberRows_];
    33183321        memset (pi,0,numberRows_*sizeof(double));
    3319         int * way = new int[numberColumns];
    3320 
     3322        int * way = new int[numberColumns_];
     3323        int numberIn = 0;
     3324       
    33213325        // Get column copy
    33223326        CoinPackedMatrix * columnCopy = matrix();
     
    33293333        const int * rowLength = copy.getVectorLengths();
    33303334        const double * elementByRow = copy.getElements();
    3331         const int * row = columnCopy->getIndices();
    3332         const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
    3333         const int * columnLength = columnCopy->getVectorLengths();
    3334         const double * element = columnCopy->getElements();
    3335 
     3335        //const int * row = columnCopy->getIndices();
     3336        //const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
     3337        //const int * columnLength = columnCopy->getVectorLengths();
     3338        //const double * element = columnCopy->getElements();
     3339       
    33363340       
    33373341        // if equality row and bounds mean artificial in basis bad
    33383342        // then do anyway
    3339 
     3343       
    33403344        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    33413345          // - if we want to reduce dj, + if we want to increase
     
    33783382          way[iColumn] = thisWay;
    33793383        }
    3380         // we need to maximize chance of doing good
    3381         int iRow;
    3382         for (iRow=0;iRow<numberRows_;iRow++) {
    3383           double lowerBound = rowLower_[iRow];
    3384           double upperBound = rowUpper_[iRow];
    3385           if (getRowStatus(iRow)==basic) {
    3386             // see if we can find a column to pivot on
    3387             int j;
    3388             double maximumDown = DBL_MAX;
    3389             double maximumUp = DBL_MAX;
    3390             int numberBad=0;
    3391             if (lowerBound<-1.0e20)
    3392               maximumDown = -1.0;;
    3393             if (upperBound<1.0e20)
    3394               maximumUp = -1.0;;
    3395             for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
    3396               int iColumn = column[j];
    3397               /* way -
    3398                  -3 - okay at upper bound with negative dj
    3399                  -2 - marginal at upper bound with zero dj - can only decrease
    3400                  -1 - bad at upper bound
    3401                  0 - we can never pivot on this row
    3402                  1 - bad at lower bound
    3403                  2 - marginal at lower bound with zero dj - can only increase
    3404                  3 - okay at lower bound with positive dj
    3405                  100 - fine we can just ignore
    3406               */
    3407               if (way[iColumn]!=100) {
    3408                 switch(way[iColumn]) {
    3409              
    3410                 case -3:
    3411                  
    3412                   break;
    3413                 case -2:
    3414                  
    3415                   break;
    3416                 case -1:
    3417                  
    3418                   break;
    3419                 case 0:
    3420                   maximumDown = -1.0;
    3421                   maximumUp=-1.0;
    3422                   break;
    3423                 case 1:
    3424                  
    3425                   break;
    3426                 case 2:
    3427                  
    3428                   break;
    3429                 case 3:
    3430                  
    3431                   break;
    3432                 default:
    3433                   break;
     3384        /*if (!numberBad)
     3385          printf("Was dual feasible before passes - rows %d\n",
     3386          numberRows_);*/
     3387        int lastNumberIn = -100000;
     3388        int numberPasses=5;
     3389        while (numberIn>lastNumberIn+numberRows_/100) {
     3390          lastNumberIn = numberIn;
     3391          // we need to maximize chance of doing good
     3392          int iRow;
     3393          for (iRow=0;iRow<numberRows_;iRow++) {
     3394            double lowerBound = rowLower_[iRow];
     3395            double upperBound = rowUpper_[iRow];
     3396            if (getRowStatus(iRow)==basic) {
     3397              // see if we can find a column to pivot on
     3398              int j;
     3399              // down is amount pi can go down
     3400              double maximumDown = DBL_MAX;
     3401              double maximumUp = DBL_MAX;
     3402              double minimumDown =0.0;
     3403              double minimumUp =0.0;
     3404              int iUp=-1;
     3405              int iDown=-1;
     3406              int iUpB=-1;
     3407              int iDownB=-1;
     3408              if (lowerBound<-1.0e20)
     3409                maximumDown = -1.0;
     3410              if (upperBound>1.0e20)
     3411                maximumUp = -1.0;
     3412              for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
     3413                int iColumn = column[j];
     3414                double value = elementByRow[j];
     3415                double djValue = dj[iColumn];
     3416                /* way -
     3417                   -3 - okay at upper bound with negative dj
     3418                   -2 - marginal at upper bound with zero dj - can only decrease
     3419                   -1 - bad at upper bound
     3420                   0 - we can never pivot on this row
     3421                   1 - bad at lower bound
     3422                   2 - marginal at lower bound with zero dj - can only increase
     3423                   3 - okay at lower bound with positive dj
     3424                   100 - fine we can just ignore
     3425                */
     3426                if (way[iColumn]!=100) {
     3427                  switch(way[iColumn]) {
     3428                   
     3429                  case -3:
     3430                    if (value>0.0) {
     3431                      if (maximumDown*value>-djValue) {
     3432                        maximumDown = -djValue/value;
     3433                        iDown = iColumn;
     3434                      }
     3435                    } else {
     3436                      if (-maximumUp*value>-djValue) {
     3437                        maximumUp = djValue/value;
     3438                        iUp = iColumn;
     3439                      }
     3440                    }
     3441                    break;
     3442                  case -2:
     3443                    if (value>0.0) {
     3444                      maximumDown = 0.0;
     3445                    } else {
     3446                      maximumUp = 0.0;
     3447                    }
     3448                    break;
     3449                  case -1:
     3450                    // see if could be satisfied
     3451                    // dj value > 0
     3452                    if (value>0.0) {
     3453                      maximumDown=0.0;
     3454                      if (maximumUp*value<djValue-dualTolerance) {
     3455                        maximumUp = 0.0; // would improve but not enough
     3456                      } else {
     3457                        if (minimumUp*value<djValue) {
     3458                          minimumUp = djValue/value;
     3459                          iUpB = iColumn;
     3460                        }
     3461                      }
     3462                    } else {
     3463                      maximumUp=0.0;
     3464                      if (-maximumDown*value<djValue-dualTolerance) {
     3465                        maximumDown = 0.0; // would improve but not enough
     3466                      } else {
     3467                        if (-minimumDown*value<djValue) {
     3468                          minimumDown = -djValue/value;
     3469                          iDownB = iColumn;
     3470                        }
     3471                      }
     3472                    }
     3473                   
     3474                    break;
     3475                  case 0:
     3476                    maximumDown = -1.0;
     3477                    maximumUp=-1.0;
     3478                    break;
     3479                  case 1:
     3480                    // see if could be satisfied
     3481                    // dj value < 0
     3482                    if (value>0.0) {
     3483                      maximumUp=0.0;
     3484                      if (maximumDown*value<-djValue-dualTolerance) {
     3485                        maximumDown = 0.0; // would improve but not enough
     3486                      } else {
     3487                        if (minimumDown*value<-djValue) {
     3488                          minimumDown = -djValue/value;
     3489                          iDownB = iColumn;
     3490                        }
     3491                      }
     3492                    } else {
     3493                      maximumDown=0.0;
     3494                      if (-maximumUp*value<-djValue-dualTolerance) {
     3495                        maximumUp = 0.0; // would improve but not enough
     3496                      } else {
     3497                        if (-minimumUp*value<-djValue) {
     3498                          minimumUp = djValue/value;
     3499                          iUpB = iColumn;
     3500                        }
     3501                      }
     3502                    }
     3503                   
     3504                    break;
     3505                  case 2:
     3506                    if (value>0.0) {
     3507                      maximumUp = 0.0;
     3508                    } else {
     3509                      maximumDown = 0.0;
     3510                    }
     3511                   
     3512                    break;
     3513                  case 3:
     3514                    if (value>0.0) {
     3515                      if (maximumUp*value>djValue) {
     3516                        maximumUp = djValue/value;
     3517                        iUp = iColumn;
     3518                      }
     3519                    } else {
     3520                      if (-maximumDown*value>djValue) {
     3521                        maximumDown = -djValue/value;
     3522                        iDown = iColumn;
     3523                      }
     3524                    }
     3525                   
     3526                    break;
     3527                  default:
     3528                    break;
     3529                  }
     3530                }
     3531              }
     3532              if (iUpB>=0)
     3533                iUp=iUpB;
     3534              if (maximumUp<=dualTolerance||maximumUp<minimumUp)
     3535                iUp=-1;
     3536              if (iDownB>=0)
     3537                iDown=iDownB;
     3538              if (maximumDown<=dualTolerance||maximumDown<minimumDown)
     3539                iDown=-1;
     3540              if (iUp>=0||iDown>=0) {
     3541                // do something
     3542                if (iUp>=0&&iDown>=0) {
     3543                  if (maximumDown>maximumUp)
     3544                    iUp=-1;
     3545                }
     3546                double change;
     3547                int kColumn;
     3548                if (iUp>=0) {
     3549                  kColumn=iUp;
     3550                  change=maximumUp;
     3551                  setRowStatus(iRow,atUpperBound);
     3552                } else {
     3553                  kColumn=iDown;
     3554                  change=-maximumDown;
     3555                  setRowStatus(iRow,atLowerBound);
     3556                }
     3557                setColumnStatus(kColumn,basic);
     3558                numberIn++;
     3559                pi[iRow]=change;
     3560                for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
     3561                  int iColumn = column[j];
     3562                  double value = elementByRow[j];
     3563                  double djValue = dj[iColumn]-change*value;
     3564                  dj[iColumn]=djValue;
     3565                  if (abs(way[iColumn])==1) {
     3566                    numberBad--;
     3567                    /*if (!numberBad)
     3568                      printf("Became dual feasible at row %d out of %d\n",
     3569                      iRow, numberRows_);*/
     3570                    lastNumberIn=-1000000;
     3571                  }
     3572                  int thisWay = 100;
     3573                  double lowerBound = columnLower_[iColumn];
     3574                  double upperBound = columnUpper_[iColumn];
     3575                  if (upperBound>lowerBound) {
     3576                    switch(getColumnStatus(iColumn)) {
     3577                     
     3578                    case basic:
     3579                      thisWay=0;
     3580                      break;
     3581                    case isFree:
     3582                    case superBasic:
     3583                      if (djValue<-dualTolerance)
     3584                        thisWay = 1;
     3585                      else if (djValue>dualTolerance)
     3586                        thisWay = -1;
     3587                      else
     3588                        { thisWay =0; abort();}
     3589                      break;
     3590                    case atUpperBound:
     3591                      if (djValue>dualTolerance)
     3592                        { thisWay =-1; abort();}
     3593                      else if (djValue<-dualTolerance)
     3594                        thisWay = -3;
     3595                      else
     3596                        thisWay = -2;
     3597                      break;
     3598                    case atLowerBound:
     3599                      if (djValue<-dualTolerance)
     3600                        { thisWay =1; abort();}
     3601                      else if (djValue>dualTolerance)
     3602                        thisWay = 3;
     3603                      else
     3604                        thisWay = 2;
     3605                      break;
     3606                    }
     3607                  }
     3608                  way[iColumn] = thisWay;
    34343609                }
    34353610              }
    34363611            }
    3437             if (max(maximumUp,maximumDown)>0.0) {
    3438               // do something
     3612          }
     3613          if (numberIn==lastNumberIn||numberBad||pivot<2)
     3614            break;
     3615          if (!(--numberPasses))
     3616            break;
     3617          //printf("%d put in so far\n",numberIn);
     3618        }
     3619        // last attempt to flip
     3620        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     3621          double lowerBound = columnLower_[iColumn];
     3622          double upperBound = columnUpper_[iColumn];
     3623          if (upperBound-lowerBound<=gap&&upperBound>lowerBound) {
     3624            double djValue=dj[iColumn];
     3625            switch(getColumnStatus(iColumn)) {
     3626             
     3627            case basic:
     3628              break;
     3629            case isFree:
     3630            case superBasic:
     3631              break;
     3632            case atUpperBound:
     3633              if (djValue>dualTolerance) {
     3634                setColumnStatus(iColumn,atUpperBound);
     3635                solution[iColumn]=upperBound;
     3636              }
     3637              break;
     3638            case atLowerBound:
     3639              if (djValue<-dualTolerance) {
     3640                setColumnStatus(iColumn,atUpperBound);
     3641                solution[iColumn]=upperBound;
     3642              }
     3643              break;
    34393644            }
    34403645          }
    34413646        }
    3442         abort();
    34433647        delete [] pi;
    34443648        delete [] dj;
    34453649        delete [] way;
    3446 #endif
     3650        handler_->message(CLP_CRASH,messages_)
     3651          <<numberIn
     3652          <<numberBad
     3653          <<CoinMessageEol;
    34473654        return -1;
    34483655      }
  • trunk/Presolve.cpp

    r59 r102  
    715715          if (!hinrow[i])
    716716            numberDropped++;
    717         printf("%d rows dropped after pass %d\n",numberDropped,
    718                iLoop+1);
     717        //printf("%d rows dropped after pass %d\n",numberDropped,
     718        //     iLoop+1);
    719719        if (numberDropped==lastDropped)
    720720          break;
  • trunk/README

    r64 r102  
    1818If you want to stress the code you can set various stuff e.g. dantzig pricing
    1919 and then go into netlib testing.  I do not guarantee that it will solve all
    20 netlib if you get too creative.
     20netlib if you get too creative.  For instance using presolve makes netlib
     21solve faster - but pilot87 prefers a large infeasibility weight.  So
    2122
    22 There will be samples in ./Samples.  Use the corresponding Makefile to
     23clp -presolve on -dualbound 1.0e10 -netlib
     24
     25works well.
     26
     27There are samples in ./Samples.  Use the corresponding Makefile to
    2328create an executable - testit.
    2429
    25 At present there are only two useful samples.
     30At present there are only three useful samples.
    2631
    2732minimum.cpp  This is the simplest possible program to read an mps file.
     33
    2834defaults.cpp.  This does not do much more, but it does it in much more
    2935complicated way by specifically setting defaults so it does give more
    3036useful information.  It also prints a solution in a format "similar" to that
    3137of MPSX.
     38
     39presolve.cpp.  This is a good driver for larger problems.
  • trunk/Test/ClpMain.cpp

    r99 r102  
    1919#include <unistd.h>
    2020#endif
    21 #define CLPVERSION "0.94"
     21#define CLPVERSION "0.94.1"
    2222
    2323//#include "CoinPackedMatrix.hpp"
     
    8181 
    8282  DIRECTION=201,DUALPIVOT,SCALING,ERRORSALLOWED,KEEPNAMES,SPARSEFACTOR,
    83   PRIMALPIVOT,PRESOLVE,
     83  PRIMALPIVOT,PRESOLVE,CRASH,
    8484 
    8585  DIRECTORY=301,IMPORT,EXPORT,RESTORE,SAVE,DUALSIMPLEX,PRIMALSIMPLEX,BAB,
     
    802802      ClpItem("presolve","Whether to presolve problem",
    803803              "off",PRESOLVE);
     804    parameters[numberParameters-1].append("on");
     805    parameters[numberParameters++]=
     806      ClpItem("crash","Whether to create basis for problem",
     807              "off",CRASH);
    804808    parameters[numberParameters-1].append("on");
    805809    parameters[numberParameters++]=
     
    10761080              preSolve = action*5;
    10771081              break;
     1082            case CRASH:
     1083              doIdiot=-1;
     1084              break;
    10781085            default:
    10791086              abort();
     
    11181125#endif
    11191126              if (type==DUALSIMPLEX) {
     1127                if (doIdiot<0)
     1128                  model2->crash(1000,2);
    11201129                model2->dual();
    11211130              } else {
    11221131#ifdef CLP_IDIOT
    1123                 if (doIdiot) {
     1132                if (doIdiot>0) {
    11241133                  Idiot info(*model2);
    11251134                  info.crash(doIdiot);
  • trunk/Test/unitTest.cpp

    r96 r102  
    301301            std::cout<<"** Analysis indicates model infeasible"
    302302                     <<std::endl;
     303          if (doIdiot<0)
     304            model2->crash(1000,2);
    303305          model2->dual();
    304306        } else {
    305307#ifdef CLP_IDIOT
    306           if (doIdiot) {
     308          if (doIdiot>0) {
    307309            Idiot info(*model2);
    308310            info.crash(doIdiot);
     
    339341#else
    340342        if (doDual) {
     343          if (doIdiot<0)
     344            solution.crash(1000,2);
    341345          solution.dual();
    342346        } else {
    343347#ifdef CLP_IDIOT
    344           if (doIdiot) {
     348          if (doIdiot>0) {
    345349            Idiot info(solution);
    346350            info.crash(doIdiot);
     
    352356      } else {
    353357        if (doDual) {
     358          if (doIdiot<0)
     359            solution.crash(1000,2);
    354360          solution.dual();
    355361        } else {
    356362#ifdef CLP_IDIOT
    357           if (doIdiot) {
     363          if (doIdiot>0) {
    358364            Idiot info(solution);
    359365            info.crash(doIdiot);
  • trunk/include/ClpMessage.hpp

    r93 r102  
    7171  CLP_IMPORT_ERRORS,
    7272  CLP_EMPTY_PROBLEM,
     73  CLP_CRASH,
    7374  CLP_DUMMY_END
    7475};
Note: See TracChangeset for help on using the changeset viewer.