Changeset 1326 for trunk


Ignore:
Timestamp:
Jan 20, 2009 2:19:35 PM (11 years ago)
Author:
forrest
Message:

First attempt at Benders decomposition

Location:
trunk/Clp
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/configure

    r1264 r1326  
    35053505ac_compiler_gnu=$ac_cv_cxx_compiler_gnu
    35063506
    3507 if test -z "$CXX" ; then
    3508   { { echo "$as_me:$LINENO: error: Failed to find a C++ compiler!" >&5
    3509 echo "$as_me: error: Failed to find a C++ compiler!" >&2;}
     3507
     3508#AC_PROG_CXX sets CXX to g++ if it cannot find a working C++ compiler
     3509#thus, we test here whether $CXX is actually working
     3510ac_ext=cc
     3511ac_cpp='$CXXCPP $CPPFLAGS'
     3512ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5'
     3513ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5'
     3514ac_compiler_gnu=$ac_cv_cxx_compiler_gnu
     3515
     3516echo "$as_me:$LINENO: checking whether C++ compiler $CXX works" >&5
     3517echo $ECHO_N "checking whether C++ compiler $CXX works... $ECHO_C" >&6;
     3518cat >conftest.$ac_ext <<_ACEOF
     3519/* confdefs.h.  */
     3520_ACEOF
     3521cat confdefs.h >>conftest.$ac_ext
     3522cat >>conftest.$ac_ext <<_ACEOF
     3523/* end confdefs.h.  */
     3524
     3525int
     3526main ()
     3527{
     3528int i=0;
     3529  ;
     3530  return 0;
     3531}
     3532_ACEOF
     3533rm -f conftest.$ac_objext
     3534if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5
     3535  (eval $ac_compile) 2>conftest.er1
     3536  ac_status=$?
     3537  grep -v '^ *+' conftest.er1 >conftest.err
     3538  rm -f conftest.er1
     3539  cat conftest.err >&5
     3540  echo "$as_me:$LINENO: \$? = $ac_status" >&5
     3541  (exit $ac_status); } &&
     3542         { ac_try='test -z "$ac_cxx_werror_flag"
     3543                         || test ! -s conftest.err'
     3544  { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5
     3545  (eval $ac_try) 2>&5
     3546  ac_status=$?
     3547  echo "$as_me:$LINENO: \$? = $ac_status" >&5
     3548  (exit $ac_status); }; } &&
     3549         { ac_try='test -s conftest.$ac_objext'
     3550  { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5
     3551  (eval $ac_try) 2>&5
     3552  ac_status=$?
     3553  echo "$as_me:$LINENO: \$? = $ac_status" >&5
     3554  (exit $ac_status); }; }; then
     3555  echo "$as_me:$LINENO: result: yes" >&5
     3556echo "${ECHO_T}yes" >&6
     3557else
     3558  echo "$as_me: failed program was:" >&5
     3559sed 's/^/| /' conftest.$ac_ext >&5
     3560
     3561echo "$as_me:$LINENO: result: no" >&5
     3562echo "${ECHO_T}no" >&6
     3563   { { echo "$as_me:$LINENO: error: failed to find a C++ compiler or C++ compiler $CXX does not work" >&5
     3564echo "$as_me: error: failed to find a C++ compiler or C++ compiler $CXX does not work" >&2;}
    35103565   { (exit 1); exit 1; }; }
    3511 fi
     3566
     3567fi
     3568rm -f conftest.err conftest.$ac_objext conftest.$ac_ext
     3569ac_ext=cc
     3570ac_cpp='$CXXCPP $CPPFLAGS'
     3571ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5'
     3572ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5'
     3573ac_compiler_gnu=$ac_cv_cxx_compiler_gnu
     3574
    35123575
    35133576# It seems that we need to cleanup something here for the Windows
     
    57165779*-*-irix6*)
    57175780  # Find out which ABI we are using.
    5718   echo '#line 5718 "configure"' > conftest.$ac_ext
     5781  echo '#line 5781 "configure"' > conftest.$ac_ext
    57195782  if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5
    57205783  (eval $ac_compile) 2>&5
     
    68506913
    68516914# Provide some information about the compiler.
    6852 echo "$as_me:6852:" \
     6915echo "$as_me:6915:" \
    68536916     "checking for Fortran 77 compiler version" >&5
    68546917ac_compiler=`set X $ac_compile; echo $2`
     
    79177980   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    79187981   -e 's:$: $lt_compiler_flag:'`
    7919    (eval echo "\"\$as_me:7919: $lt_compile\"" >&5)
     7982   (eval echo "\"\$as_me:7982: $lt_compile\"" >&5)
    79207983   (eval "$lt_compile" 2>conftest.err)
    79217984   ac_status=$?
    79227985   cat conftest.err >&5
    7923    echo "$as_me:7923: \$? = $ac_status" >&5
     7986   echo "$as_me:7986: \$? = $ac_status" >&5
    79247987   if (exit $ac_status) && test -s "$ac_outfile"; then
    79257988     # The compiler can only warn and ignore the option if not recognized
     
    81858248   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    81868249   -e 's:$: $lt_compiler_flag:'`
    8187    (eval echo "\"\$as_me:8187: $lt_compile\"" >&5)
     8250   (eval echo "\"\$as_me:8250: $lt_compile\"" >&5)
    81888251   (eval "$lt_compile" 2>conftest.err)
    81898252   ac_status=$?
    81908253   cat conftest.err >&5
    8191    echo "$as_me:8191: \$? = $ac_status" >&5
     8254   echo "$as_me:8254: \$? = $ac_status" >&5
    81928255   if (exit $ac_status) && test -s "$ac_outfile"; then
    81938256     # The compiler can only warn and ignore the option if not recognized
     
    82898352   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    82908353   -e 's:$: $lt_compiler_flag:'`
    8291    (eval echo "\"\$as_me:8291: $lt_compile\"" >&5)
     8354   (eval echo "\"\$as_me:8354: $lt_compile\"" >&5)
    82928355   (eval "$lt_compile" 2>out/conftest.err)
    82938356   ac_status=$?
    82948357   cat out/conftest.err >&5
    8295    echo "$as_me:8295: \$? = $ac_status" >&5
     8358   echo "$as_me:8358: \$? = $ac_status" >&5
    82968359   if (exit $ac_status) && test -s out/conftest2.$ac_objext
    82978360   then
     
    1063410697  lt_status=$lt_dlunknown
    1063510698  cat > conftest.$ac_ext <<EOF
    10636 #line 10636 "configure"
     10699#line 10699 "configure"
    1063710700#include "confdefs.h"
    1063810701
     
    1073410797  lt_status=$lt_dlunknown
    1073510798  cat > conftest.$ac_ext <<EOF
    10736 #line 10736 "configure"
     10799#line 10799 "configure"
    1073710800#include "confdefs.h"
    1073810801
     
    1307813141   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    1307913142   -e 's:$: $lt_compiler_flag:'`
    13080    (eval echo "\"\$as_me:13080: $lt_compile\"" >&5)
     13143   (eval echo "\"\$as_me:13143: $lt_compile\"" >&5)
    1308113144   (eval "$lt_compile" 2>conftest.err)
    1308213145   ac_status=$?
    1308313146   cat conftest.err >&5
    13084    echo "$as_me:13084: \$? = $ac_status" >&5
     13147   echo "$as_me:13147: \$? = $ac_status" >&5
    1308513148   if (exit $ac_status) && test -s "$ac_outfile"; then
    1308613149     # The compiler can only warn and ignore the option if not recognized
     
    1318213245   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    1318313246   -e 's:$: $lt_compiler_flag:'`
    13184    (eval echo "\"\$as_me:13184: $lt_compile\"" >&5)
     13247   (eval echo "\"\$as_me:13247: $lt_compile\"" >&5)
    1318513248   (eval "$lt_compile" 2>out/conftest.err)
    1318613249   ac_status=$?
    1318713250   cat out/conftest.err >&5
    13188    echo "$as_me:13188: \$? = $ac_status" >&5
     13251   echo "$as_me:13251: \$? = $ac_status" >&5
    1318913252   if (exit $ac_status) && test -s out/conftest2.$ac_objext
    1319013253   then
     
    1475214815   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    1475314816   -e 's:$: $lt_compiler_flag:'`
    14754    (eval echo "\"\$as_me:14754: $lt_compile\"" >&5)
     14817   (eval echo "\"\$as_me:14817: $lt_compile\"" >&5)
    1475514818   (eval "$lt_compile" 2>conftest.err)
    1475614819   ac_status=$?
    1475714820   cat conftest.err >&5
    14758    echo "$as_me:14758: \$? = $ac_status" >&5
     14821   echo "$as_me:14821: \$? = $ac_status" >&5
    1475914822   if (exit $ac_status) && test -s "$ac_outfile"; then
    1476014823     # The compiler can only warn and ignore the option if not recognized
     
    1485614919   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    1485714920   -e 's:$: $lt_compiler_flag:'`
    14858    (eval echo "\"\$as_me:14858: $lt_compile\"" >&5)
     14921   (eval echo "\"\$as_me:14921: $lt_compile\"" >&5)
    1485914922   (eval "$lt_compile" 2>out/conftest.err)
    1486014923   ac_status=$?
    1486114924   cat out/conftest.err >&5
    14862    echo "$as_me:14862: \$? = $ac_status" >&5
     14925   echo "$as_me:14925: \$? = $ac_status" >&5
    1486314926   if (exit $ac_status) && test -s out/conftest2.$ac_objext
    1486414927   then
     
    1706317126   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    1706417127   -e 's:$: $lt_compiler_flag:'`
    17065    (eval echo "\"\$as_me:17065: $lt_compile\"" >&5)
     17128   (eval echo "\"\$as_me:17128: $lt_compile\"" >&5)
    1706617129   (eval "$lt_compile" 2>conftest.err)
    1706717130   ac_status=$?
    1706817131   cat conftest.err >&5
    17069    echo "$as_me:17069: \$? = $ac_status" >&5
     17132   echo "$as_me:17132: \$? = $ac_status" >&5
    1707017133   if (exit $ac_status) && test -s "$ac_outfile"; then
    1707117134     # The compiler can only warn and ignore the option if not recognized
     
    1733117394   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    1733217395   -e 's:$: $lt_compiler_flag:'`
    17333    (eval echo "\"\$as_me:17333: $lt_compile\"" >&5)
     17396   (eval echo "\"\$as_me:17396: $lt_compile\"" >&5)
    1733417397   (eval "$lt_compile" 2>conftest.err)
    1733517398   ac_status=$?
    1733617399   cat conftest.err >&5
    17337    echo "$as_me:17337: \$? = $ac_status" >&5
     17400   echo "$as_me:17400: \$? = $ac_status" >&5
    1733817401   if (exit $ac_status) && test -s "$ac_outfile"; then
    1733917402     # The compiler can only warn and ignore the option if not recognized
     
    1743517498   -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \
    1743617499   -e 's:$: $lt_compiler_flag:'`
    17437    (eval echo "\"\$as_me:17437: $lt_compile\"" >&5)
     17500   (eval echo "\"\$as_me:17500: $lt_compile\"" >&5)
    1743817501   (eval "$lt_compile" 2>out/conftest.err)
    1743917502   ac_status=$?
    1744017503   cat out/conftest.err >&5
    17441    echo "$as_me:17441: \$? = $ac_status" >&5
     17504   echo "$as_me:17504: \$? = $ac_status" >&5
    1744217505   if (exit $ac_status) && test -s out/conftest2.$ac_objext
    1744317506   then
  • trunk/Clp/src/ClpSimplex.hpp

    r1321 r1326  
    715715  */
    716716  double scaleObjective(double value);
     717  /// Solve using Dantzig-Wolfe decomposition and maybe in parallel
     718  int solveDW(CoinStructuredModel * model);
     719  /// Solve using Benders decomposition and maybe in parallel
     720  int solveBenders(CoinStructuredModel * model);
    717721public:
    718722  /** For advanced use.  When doing iterative solves things can get
  • trunk/Clp/src/ClpSolve.cpp

    r1325 r1326  
    31253125ClpSimplex::solve(CoinStructuredModel * model)
    31263126{
    3127   double time1 = CoinCpuTime();
    31283127  // analyze structure
    3129   int numberColumns = model->numberColumns();
    31303128  int numberRowBlocks=model->numberRowBlocks();
    31313129  int numberColumnBlocks = model->numberColumnBlocks();
     
    31643162    }
    31653163  }
    3166   if (numberRowBlocks==numberColumnBlocks||numberRowBlocks==numberColumnBlocks+1) {
    3167     // looks like Dantzig-Wolfe
     3164  int * rowCounts = new int [numberRowBlocks];
     3165  CoinZeroN(rowCounts,numberRowBlocks);
     3166  int * columnCounts = new int [numberColumnBlocks+1];
     3167  CoinZeroN(columnCounts,numberColumnBlocks);
     3168  int decomposeType=0;
     3169  for (int i=0;i<numberElementBlocks;i++) {
     3170    int iRowBlock = blockInfo[i].rowBlock;
     3171    int iColumnBlock = blockInfo[i].columnBlock;
     3172    rowCounts[iRowBlock]++;
     3173    columnCounts[iColumnBlock]++;
     3174  }
     3175  if (numberRowBlocks==numberColumnBlocks||
     3176      numberRowBlocks==numberColumnBlocks+1) {
     3177    // could be Dantzig-Wolfe
     3178    int numberG1=0;
     3179    for (int i=0;i<numberRowBlocks;i++) {
     3180      if (rowCounts[i]>1)
     3181        numberG1++;
     3182    }
    31683183    bool masterColumns = (numberColumnBlocks==numberRowBlocks);
    31693184    if ((masterColumns&&numberElementBlocks==2*numberRowBlocks-1)
    31703185        ||(!masterColumns&&numberElementBlocks==2*numberRowBlocks)) {
    3171       // make up problems
    3172       int numberBlocks=numberRowBlocks-1;
    3173       // Find master rows and columns
    3174       int * rowCounts = new int [numberRowBlocks];
    3175       CoinZeroN(rowCounts,numberRowBlocks);
    3176       int * columnCounts = new int [numberColumnBlocks+1];
    3177       CoinZeroN(columnCounts,numberColumnBlocks);
    3178       int iBlock;
    3179       for (iBlock = 0;iBlock<numberElementBlocks;iBlock++) {
    3180         int iRowBlock = blockInfo[iBlock].rowBlock;
    3181         rowCounts[iRowBlock]++;
    3182         int iColumnBlock =blockInfo[iBlock].columnBlock;
    3183         columnCounts[iColumnBlock]++;
    3184       }
    3185       int * whichBlock = new int [numberElementBlocks];
    3186       int masterRowBlock=-1;
    3187       for (iBlock = 0;iBlock<numberRowBlocks;iBlock++) {
    3188         if (rowCounts[iBlock]>1) {
    3189           if (masterRowBlock==-1) {
    3190             masterRowBlock=iBlock;
    3191           } else {
    3192             // Can't decode
    3193             masterRowBlock=-2;
    3194             break;
    3195           }
    3196         }
    3197       }
    3198       int masterColumnBlock=-1;
    3199       int kBlock=0;
    3200       for (iBlock = 0;iBlock<numberColumnBlocks;iBlock++) {
    3201         int count=columnCounts[iBlock];
    3202         columnCounts[iBlock]=kBlock;
    3203         kBlock += count;
    3204       }
    3205       for (iBlock = 0;iBlock<numberElementBlocks;iBlock++) {
    3206         int iColumnBlock = blockInfo[iBlock].columnBlock;
    3207         whichBlock[columnCounts[iColumnBlock]]=iBlock;
    3208         columnCounts[iColumnBlock]++;
    3209       }
    3210       for (iBlock = numberColumnBlocks-1;iBlock>=0;iBlock--)
    3211         columnCounts[iBlock+1]=columnCounts[iBlock];
    3212       columnCounts[0]=0;
    3213       for (iBlock = 0;iBlock<numberColumnBlocks;iBlock++) {
    3214         int count=columnCounts[iBlock+1]-columnCounts[iBlock];
    3215         if (count==1) {
    3216           int kBlock = whichBlock[columnCounts[iBlock]];
    3217           int iRowBlock = blockInfo[kBlock].rowBlock;
    3218           if (iRowBlock==masterRowBlock) {
    3219             if (masterColumnBlock==-1) {
    3220               masterColumnBlock=iBlock;
    3221             } else {
    3222               // Can't decode
    3223               masterColumnBlock=-2;
    3224               break;
    3225             }
    3226           }
    3227         }
    3228       }
    3229       if (masterRowBlock<0||masterColumnBlock==-2) {
    3230         // What now
    3231         abort();
    3232       }
    3233       delete [] rowCounts;
    3234       // create all data
    3235       const CoinPackedMatrix ** top = new const CoinPackedMatrix * [numberColumnBlocks];
    3236       ClpSimplex * sub = new ClpSimplex [numberBlocks];
    3237       ClpSimplex master;
    3238       // Set offset
    3239       master.setObjectiveOffset(model->objectiveOffset());
    3240       kBlock=0;
    3241       int masterBlock=-1;
    3242       for (iBlock = 0;iBlock<numberColumnBlocks;iBlock++) {
    3243         top[kBlock]=NULL;
    3244         int start=columnCounts[iBlock];
    3245         int end=columnCounts[iBlock+1];
    3246         assert (end-start<=2);
    3247         for (int j=start;j<end;j++) {
    3248           int jBlock = whichBlock[j];
    3249           int iRowBlock = blockInfo[jBlock].rowBlock;
    3250           int iColumnBlock =blockInfo[jBlock].columnBlock;
    3251           assert (iColumnBlock==iBlock);
    3252           if (iColumnBlock!=masterColumnBlock&&iRowBlock==masterRowBlock) {
    3253             // top matrix
    3254             top[kBlock]=model->coinBlock(jBlock)->packedMatrix();
    3255           } else {
    3256             const CoinPackedMatrix * matrix
    3257               = model->coinBlock(jBlock)->packedMatrix();
    3258             // Get pointers to arrays
    3259             const double * rowLower;
    3260             const double * rowUpper;
    3261             const double * columnLower;
    3262             const double * columnUpper;
    3263             const double * objective;
    3264             model->block(iRowBlock,iColumnBlock,rowLower,rowUpper,
    3265                          columnLower,columnUpper,objective);
    3266             if (iColumnBlock!=masterColumnBlock) {
    3267               // diagonal block
    3268               sub[kBlock].loadProblem(*matrix,columnLower,columnUpper,
    3269                                 objective,rowLower,rowUpper);
    3270               if (true) {
    3271                 double * lower = sub[kBlock].columnLower();
    3272                 double * upper = sub[kBlock].columnUpper();
    3273                 int n = sub[kBlock].numberColumns();
    3274                 for (int i=0;i<n;i++) {
    3275                   lower[i]=CoinMax(-1.0e8,lower[i]);
    3276                   upper[i]=CoinMin(1.0e8,upper[i]);
    3277                 }
    3278               }
    3279               if (optimizationDirection_<0.0) {
    3280                 double * obj = sub[kBlock].objective();
    3281                 int n = sub[kBlock].numberColumns();
    3282                 for (int i=0;i<n;i++)
    3283                   obj[i] = - obj[i];
    3284               }
    3285               if (this->factorizationFrequency()==200) {
    3286                 // User did not touch preset
    3287                 sub[kBlock].defaultFactorizationFrequency();
    3288               } else {
    3289                 // make sure model has correct value
    3290                 sub[kBlock].setFactorizationFrequency(this->factorizationFrequency());
    3291               }
    3292               sub[kBlock].setPerturbation(50);
    3293               // Set columnCounts to be diagonal block index for cleanup
    3294               columnCounts[kBlock]=jBlock;
    3295             } else {
    3296               // master
    3297               masterBlock=jBlock;
    3298               master.loadProblem(*matrix,columnLower,columnUpper,
    3299                                 objective,rowLower,rowUpper);
    3300               if (optimizationDirection_<0.0) {
    3301                 double * obj = master.objective();
    3302                 int n = master.numberColumns();
    3303                 for (int i=0;i<n;i++)
    3304                   obj[i] = - obj[i];
    3305               }
    3306             }
    3307           }
    3308         }
    3309         if (iBlock!=masterColumnBlock)
    3310           kBlock++;
    3311       }
    3312       delete [] whichBlock;
    3313       delete [] blockInfo;
    3314       // For now master must have been defined (does not have to have columns)
    3315       assert (master.numberRows());
    3316       assert (masterBlock>=0);
    3317       int numberMasterRows = master.numberRows();
    3318       // Overkill in terms of space
    3319       int spaceNeeded = CoinMax(numberBlocks*(numberMasterRows+1),
    3320                                 2*numberMasterRows);
    3321       int * rowAdd = new int[spaceNeeded];
    3322       double * elementAdd = new double[spaceNeeded];
    3323       spaceNeeded = numberBlocks;
    3324       int * columnAdd = new int[spaceNeeded+1];
    3325       double * objective = new double[spaceNeeded];
    3326       // Add in costed slacks
    3327       int firstArtificial = master.numberColumns();
    3328       int lastArtificial = firstArtificial;
    3329       if (true) {
    3330         const double * lower = master.rowLower();
    3331         const double * upper = master.rowUpper();
    3332         int kCol=0;
    3333         for (int iRow=0;iRow<numberMasterRows;iRow++) {
    3334           if (lower[iRow]>-1.0e10) {
    3335             rowAdd[kCol]=iRow;
    3336             elementAdd[kCol++]=1.0;
    3337           }
    3338           if (upper[iRow]<1.0e10) {
    3339             rowAdd[kCol]=iRow;
    3340             elementAdd[kCol++]=-1.0;
    3341           }
    3342         }
    3343         if (kCol>spaceNeeded) {
    3344           spaceNeeded = kCol;
    3345           delete [] columnAdd;
    3346           delete [] objective;
    3347           columnAdd = new int[spaceNeeded+1];
    3348           objective = new double[spaceNeeded];
    3349         }
    3350         for (int i=0;i<kCol;i++) {
    3351           columnAdd[i]=i;
    3352           objective[i]=1.0e13;
    3353         }
    3354         columnAdd[kCol]=kCol;
    3355         master.addColumns(kCol,NULL,NULL,objective,
    3356                           columnAdd,rowAdd,elementAdd);
    3357         lastArtificial = master.numberColumns();
    3358       }
    3359       int maxPass=500;
    3360       int iPass;
    3361       double lastObjective=1.0e31;
    3362       // Create convexity rows for proposals
    3363       int numberMasterColumns = master.numberColumns();
    3364       master.resize(numberMasterRows+numberBlocks,numberMasterColumns);
    3365       if (this->factorizationFrequency()==200) {
    3366         // User did not touch preset
    3367         master.defaultFactorizationFrequency();
    3368       } else {
    3369         // make sure model has correct value
    3370         master.setFactorizationFrequency(this->factorizationFrequency());
    3371       }
    3372       master.setPerturbation(50);
    3373       // Arrays to say which block and when created
    3374       int maximumColumns = 2*numberMasterRows+10*numberBlocks;
    3375       whichBlock = new int[maximumColumns];
    3376       int * when = new int[maximumColumns];
    3377       int numberColumnsGenerated=numberBlocks;
    3378       // fill in rhs and add in artificials
    3379       {
    3380         double * rowLower = master.rowLower();
    3381         double * rowUpper = master.rowUpper();
    3382         int iBlock;
    3383         columnAdd[0] = 0;
    3384         for (iBlock=0;iBlock<numberBlocks;iBlock++) {
    3385           int iRow = iBlock + numberMasterRows;;
    3386           rowLower[iRow]=1.0;
    3387           rowUpper[iRow]=1.0;
    3388           rowAdd[iBlock] = iRow;
    3389           elementAdd[iBlock] = 1.0;
    3390           objective[iBlock] = 1.0e13;
    3391           columnAdd[iBlock+1] = iBlock+1;
    3392           when[iBlock]=-1;
    3393           whichBlock[iBlock] = iBlock;
    3394         }
    3395         master.addColumns(numberBlocks,NULL,NULL,objective,
    3396                           columnAdd, rowAdd, elementAdd);
    3397       }
    3398       // and resize matrix to double check clp will be happy
    3399       //master.matrix()->setDimensions(numberMasterRows+numberBlocks,
    3400       //                         numberMasterColumns+numberBlocks);
    3401       std::cout<<"Time to decompose "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
    3402       for (iPass=0;iPass<maxPass;iPass++) {
    3403         printf("Start of pass %d\n",iPass);
    3404         // Solve master - may be infeasible
    3405         //master.scaling(0);
    3406         if (0) {
    3407           master.writeMps("yy.mps");
    3408         }
    3409         // Correct artificials
    3410         double sumArtificials=0.0;
    3411         if (iPass) {
    3412           double * upper = master.columnUpper();
    3413           double * solution = master.primalColumnSolution();
    3414           double * obj = master.objective();
    3415           sumArtificials=0.0;
    3416           for (int i=firstArtificial;i<lastArtificial;i++) {
    3417             sumArtificials += solution[i];
    3418             //assert (solution[i]>-1.0e-2);
    3419             if (solution[i]<1.0e-6) {
    3420 #if 0
    3421               // Could take out
    3422               obj[i]=0.0;
    3423               upper[i]=0.0;
    3424 #else
    3425               obj[i]=1.0e7;
    3426               upper[i]=1.0e-1;
    3427 #endif
    3428               solution[i]=0.0;
    3429               master.setColumnStatus(i,isFixed);
    3430             } else {
    3431               upper[i]=solution[i]+1.0e-5*(1.0+solution[i]);
    3432             }
    3433           }
    3434           printf("Sum of artificials before solve is %g\n",sumArtificials);
    3435         }
    3436         // scale objective to be reasonable
    3437         double scaleFactor = master.scaleObjective(-1.0e9);
    3438         {
    3439           double * dual = master.dualRowSolution();
    3440           int n=master.numberRows();
    3441           memset(dual,0,n*sizeof(double));
    3442           double * solution = master.primalColumnSolution();
    3443           master.clpMatrix()->times(1.0,solution,dual);
    3444           double sum=0.0;
    3445           double * lower = master.rowLower();
    3446           double * upper = master.rowUpper();
    3447           for (int iRow=0;iRow<n;iRow++) {
    3448             double value = dual[iRow];
    3449             if (value>upper[iRow])
    3450               sum += value-upper[iRow];
    3451             else if (value<lower[iRow])
    3452               sum -= value-lower[iRow];
    3453           }
    3454           printf("suminf %g\n",sum);
    3455           lower = master.columnLower();
    3456           upper = master.columnUpper();
    3457           n=master.numberColumns();
    3458           for (int iColumn=0;iColumn<n;iColumn++) {
    3459             double value = solution[iColumn];
    3460             if (value>upper[iColumn]+1.0e-5)
    3461               sum += value-upper[iColumn];
    3462             else if (value<lower[iColumn]-1.0e-5)
    3463               sum -= value-lower[iColumn];
    3464           }
    3465           printf("suminf %g\n",sum);
    3466         }
    3467         master.primal(1);
    3468         // Correct artificials
    3469         sumArtificials=0.0;
    3470         {
    3471           double * solution = master.primalColumnSolution();
    3472           for (int i=firstArtificial;i<lastArtificial;i++) {
    3473             sumArtificials += solution[i];
    3474           }
    3475           printf("Sum of artificials after solve is %g\n",sumArtificials);
    3476         }
    3477         master.scaleObjective(scaleFactor);
    3478         int problemStatus = master.status(); // do here as can change (delcols)
    3479         if (master.numberIterations()==0&&iPass)
    3480           break; // finished
    3481         if (master.objectiveValue()>lastObjective-1.0e-7&&iPass>555)
    3482           break; // finished
    3483         lastObjective = master.objectiveValue();
    3484         // mark basic ones and delete if necessary
    3485         int iColumn;
    3486         numberColumnsGenerated=master.numberColumns()-numberMasterColumns;
    3487         for (iColumn=0;iColumn<numberColumnsGenerated;iColumn++) {
    3488           if (master.getStatus(iColumn+numberMasterColumns)==ClpSimplex::basic)
    3489             when[iColumn]=iPass;
    3490         }
    3491         if (numberColumnsGenerated+numberBlocks>maximumColumns) {
    3492           // delete
    3493           int numberKeep=0;
    3494           int numberDelete=0;
    3495           int * whichDelete = new int[numberColumnsGenerated];
    3496           for (iColumn=0;iColumn<numberColumnsGenerated;iColumn++) {
    3497             if (when[iColumn]>iPass-7) {
    3498               // keep
    3499               when[numberKeep]=when[iColumn];
    3500               whichBlock[numberKeep++]=whichBlock[iColumn];
    3501             } else {
    3502               // delete
    3503               whichDelete[numberDelete++]=iColumn+numberMasterColumns;
    3504             }
    3505           }
    3506           numberColumnsGenerated -= numberDelete;
    3507           master.deleteColumns(numberDelete,whichDelete);
    3508           delete [] whichDelete;
    3509         }
    3510         const double * dual=NULL;
    3511         bool deleteDual=false;
    3512         if (problemStatus==0) {
    3513           dual = master.dualRowSolution();
    3514         } else if (problemStatus==1) {
    3515           // could do composite objective
    3516           dual = master.infeasibilityRay();
    3517           deleteDual = true;
    3518           printf("The sum of infeasibilities is %g\n",
    3519                  master.sumPrimalInfeasibilities());
    3520         } else if (!master.numberColumns()) {
    3521           assert(!iPass);
    3522           dual = master.dualRowSolution();
    3523           memset(master.dualRowSolution(),
    3524                  0,(numberMasterRows+numberBlocks)*sizeof(double));
    3525         } else {
    3526           abort();
    3527         }
    3528         // Scale back on first time
    3529         if (!iPass) {
    3530           double * dual2 = master.dualRowSolution();
    3531           for (int i=0;i<numberMasterRows+numberBlocks;i++) {
    3532             dual2[i] *= 1.0e-7;
    3533           }
    3534           dual = master.dualRowSolution();
    3535         }
    3536         // Create objective for sub problems and solve
    3537         columnAdd[0]=0;
    3538         int numberProposals=0;
    3539         for (iBlock=0;iBlock<numberBlocks;iBlock++) {
    3540           int numberColumns2 = sub[iBlock].numberColumns();
    3541           double * saveObj = new double [numberColumns2];
    3542           double * objective2 = sub[iBlock].objective();
    3543           memcpy(saveObj,objective2,numberColumns2*sizeof(double));
    3544           // new objective
    3545           top[iBlock]->transposeTimes(dual,objective2);
    3546           int i;
    3547           if (problemStatus==0) {
    3548             for (i=0;i<numberColumns2;i++)
    3549               objective2[i] = saveObj[i]-objective2[i];
    3550           } else {
    3551             for (i=0;i<numberColumns2;i++)
    3552               objective2[i] = -objective2[i];
    3553           }
    3554           // scale objective to be reasonable
    3555           double scaleFactor =
    3556             sub[iBlock].scaleObjective((sumArtificials>1.0e-5) ? -1.0e-4 :-1.0e9);
    3557           if (iPass) {
    3558             sub[iBlock].primal();
    3559           } else {
    3560             sub[iBlock].dual();
    3561           }
    3562           sub[iBlock].scaleObjective(scaleFactor);
    3563           if (!sub[iBlock].isProvenOptimal()&&
    3564               !sub[iBlock].isProvenDualInfeasible()) {
    3565             memset(objective2,0,numberColumns2*sizeof(double));
    3566             sub[iBlock].primal();
    3567             if (problemStatus==0) {
    3568               for (i=0;i<numberColumns2;i++)
    3569                 objective2[i] = saveObj[i]-objective2[i];
    3570             } else {
    3571               for (i=0;i<numberColumns2;i++)
    3572                 objective2[i] = -objective2[i];
    3573             }
    3574             double scaleFactor = sub[iBlock].scaleObjective(-1.0e9);
    3575             sub[iBlock].primal(1);
    3576             sub[iBlock].scaleObjective(scaleFactor);
    3577           }
    3578           memcpy(objective2,saveObj,numberColumns2*sizeof(double));
    3579           // get proposal
    3580           if (sub[iBlock].numberIterations()||!iPass) {
    3581             double objValue =0.0;
    3582             int start = columnAdd[numberProposals];
    3583             // proposal
    3584             if (sub[iBlock].isProvenOptimal()) {
    3585               const double * solution = sub[iBlock].primalColumnSolution();
    3586               top[iBlock]->times(solution,elementAdd+start);
    3587               for (i=0;i<numberColumns2;i++)
    3588                 objValue += solution[i]*saveObj[i];
    3589               // See if good dj and pack down
    3590               int number=start;
    3591               double dj = objValue;
    3592               if (problemStatus)
    3593                 dj=0.0;
    3594               double smallest=1.0e100;
    3595               double largest=0.0;
    3596               for (i=0;i<numberMasterRows;i++) {
    3597                 double value = elementAdd[start+i];
    3598                 if (fabs(value)>1.0e-15) {
    3599                   dj -= dual[i]*value;
    3600                   smallest = CoinMin(smallest,fabs(value));
    3601                   largest = CoinMax(largest,fabs(value));
    3602                   rowAdd[number]=i;
    3603                   elementAdd[number++]=value;
    3604                 }
    3605               }
    3606               // and convexity
    3607               dj -= dual[numberMasterRows+iBlock];
    3608               rowAdd[number]=numberMasterRows+iBlock;
    3609               elementAdd[number++]=1.0;
    3610               // if elements large then scale?
    3611               //if (largest>1.0e8||smallest<1.0e-8)
    3612               printf("For subproblem %d smallest - %g, largest %g - dj %g\n",
    3613                      iBlock,smallest,largest,dj);
    3614               if (dj<-1.0e-6||!iPass) {
    3615                 // take
    3616                 objective[numberProposals]=objValue;
    3617                 columnAdd[++numberProposals]=number;
    3618                 when[numberColumnsGenerated]=iPass;
    3619                 whichBlock[numberColumnsGenerated++]=iBlock;
    3620               }
    3621             } else if (sub[iBlock].isProvenDualInfeasible()) {
    3622               // use ray
    3623               const double * solution = sub[iBlock].unboundedRay();
    3624               top[iBlock]->times(solution,elementAdd+start);
    3625               for (i=0;i<numberColumns2;i++)
    3626                 objValue += solution[i]*saveObj[i];
    3627               // See if good dj and pack down
    3628               int number=start;
    3629               double dj = objValue;
    3630               double smallest=1.0e100;
    3631               double largest=0.0;
    3632               for (i=0;i<numberMasterRows;i++) {
    3633                 double value = elementAdd[start+i];
    3634                 if (fabs(value)>1.0e-15) {
    3635                   dj -= dual[i]*value;
    3636                   smallest = CoinMin(smallest,fabs(value));
    3637                   largest = CoinMax(largest,fabs(value));
    3638                   rowAdd[number]=i;
    3639                   elementAdd[number++]=value;
    3640                 }
    3641               }
    3642               // if elements large or small then scale?
    3643               //if (largest>1.0e8||smallest<1.0e-8)
    3644               printf("For subproblem ray %d smallest - %g, largest %g - dj %g\n",
    3645                      iBlock,smallest,largest,dj);
    3646               if (dj<-1.0e-6) {
    3647                 // take
    3648                 objective[numberProposals]=objValue;
    3649                 columnAdd[++numberProposals]=number;
    3650                 when[numberColumnsGenerated]=iPass;
    3651                 whichBlock[numberColumnsGenerated++]=iBlock;
    3652               }
    3653             } else {
    3654               abort();
    3655             }
    3656           }
    3657           delete [] saveObj;
    3658         }
    3659         if (deleteDual)
    3660           delete [] dual;
    3661         if (numberProposals)
    3662           master.addColumns(numberProposals,NULL,NULL,objective,
    3663                             columnAdd,rowAdd,elementAdd);
    3664       }
    3665       std::cout<<"Time at end of D-W "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
    3666       //master.scaling(0);
    3667       //master.primal(1);
    3668       loadProblem(*model);
    3669       // now put back a good solution
    3670       double * lower = new double[numberMasterRows+numberBlocks];
    3671       double * upper = new double[numberMasterRows+numberBlocks];
    3672       numberColumnsGenerated  += numberMasterColumns;
    3673       double * sol = new double[numberColumnsGenerated];
    3674       const double * solution = master.primalColumnSolution();
    3675       const double * masterLower = master.rowLower();
    3676       const double * masterUpper = master.rowUpper();
    3677       double * fullSolution = primalColumnSolution();
    3678       const double * fullLower = columnLower();
    3679       const double * fullUpper = columnUpper();
    3680       const double * rowSolution = master.primalRowSolution();
    3681       double * fullRowSolution = primalRowSolution();
    3682       const int * rowBack = model->coinBlock(masterBlock)->originalRows();
    3683       int numberRows2 = model->coinBlock(masterBlock)->numberRows();
    3684       const int * columnBack = model->coinBlock(masterBlock)->originalColumns();
    3685       int numberColumns2 = model->coinBlock(masterBlock)->numberColumns();
    3686       for (int iRow=0;iRow<numberRows2;iRow++) {
    3687         int kRow = rowBack[iRow];
    3688         setRowStatus(kRow,master.getRowStatus(iRow));
    3689         fullRowSolution[kRow]=rowSolution[iRow];
    3690       }
    3691       for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
    3692         int kColumn = columnBack[iColumn];
    3693         setStatus(kColumn,master.getStatus(iColumn));
    3694         fullSolution[kColumn]=solution[iColumn];
    3695       }
    3696       for (iBlock=0;iBlock<numberBlocks;iBlock++) {
    3697         // move basis
    3698         int kBlock = columnCounts[iBlock];
    3699         const int * rowBack = model->coinBlock(kBlock)->originalRows();
    3700         int numberRows2 = model->coinBlock(kBlock)->numberRows();
    3701         const int * columnBack = model->coinBlock(kBlock)->originalColumns();
    3702         int numberColumns2 = model->coinBlock(kBlock)->numberColumns();
    3703         for (int iRow=0;iRow<numberRows2;iRow++) {
    3704           int kRow = rowBack[iRow];
    3705           setRowStatus(kRow,sub[iBlock].getRowStatus(iRow));
    3706         }
    3707         for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
    3708           int kColumn = columnBack[iColumn];
    3709           setStatus(kColumn,sub[iBlock].getStatus(iColumn));
    3710         }
    3711         // convert top bit to by rows
    3712         CoinPackedMatrix topMatrix = *top[iBlock];
    3713         topMatrix.reverseOrdering();
    3714         // zero solution
    3715         memset(sol,0,numberColumnsGenerated*sizeof(double));
    3716        
    3717         for (int i=numberMasterColumns;i<numberColumnsGenerated;i++) {
    3718           if (whichBlock[i-numberMasterColumns]==iBlock)
    3719             sol[i] = solution[i];
    3720         }
    3721         memset(lower,0,(numberMasterRows+numberBlocks)*sizeof(double));
    3722         master.clpMatrix()->times(1.0,sol,lower);
    3723         for (int iRow=0;iRow<numberMasterRows;iRow++) {
    3724           double value = lower[iRow];
    3725           if (masterUpper[iRow]<1.0e20)
    3726             upper[iRow] = value;
    3727           else
    3728             upper[iRow]=COIN_DBL_MAX;
    3729           if (masterLower[iRow]>-1.0e20)
    3730             lower[iRow] = value;
    3731           else
    3732             lower[iRow]=-COIN_DBL_MAX;
    3733         }
    3734         sub[iBlock].addRows(numberMasterRows,lower,upper,
    3735                             topMatrix.getVectorStarts(),
    3736                             topMatrix.getVectorLengths(),
    3737                             topMatrix.getIndices(),
    3738                             topMatrix.getElements());
    3739         sub[iBlock].primal(1);
    3740         const double * subSolution = sub[iBlock].primalColumnSolution();
    3741         const double * subRowSolution = sub[iBlock].primalRowSolution();
    3742         // move solution
    3743         for (int iRow=0;iRow<numberRows2;iRow++) {
    3744           int kRow = rowBack[iRow];
    3745           fullRowSolution[kRow]=subRowSolution[iRow];
    3746         }
    3747         for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
    3748           int kColumn = columnBack[iColumn];
    3749           fullSolution[kColumn]=subSolution[iColumn];
    3750         }
    3751       }
    3752       for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    3753         if (fullSolution[iColumn]<fullUpper[iColumn]-1.0e-8&&
    3754             fullSolution[iColumn]>fullLower[iColumn]+1.0e-8) {
    3755           if (getStatus(iColumn)!=ClpSimplex::basic) {
    3756             if (columnLower_[iColumn]>-1.0e30||
    3757                 columnUpper_[iColumn]<1.0e30)
    3758               setStatus(iColumn,ClpSimplex::superBasic);
    3759             else
    3760               setStatus(iColumn,ClpSimplex::isFree);
    3761           }
    3762         } else if (fullSolution[iColumn]>=fullUpper[iColumn]-1.0e-8) {
    3763           // may help to make rest non basic
    3764           if (getStatus(iColumn)!=ClpSimplex::basic)
    3765             setStatus(iColumn,ClpSimplex::atUpperBound);
    3766         } else if (fullSolution[iColumn]<=fullLower[iColumn]+1.0e-8) {
    3767           // may help to make rest non basic
    3768           if (getStatus(iColumn)!=ClpSimplex::basic)
    3769             setStatus(iColumn,ClpSimplex::atLowerBound);
    3770         }
    3771       }
    3772       //int numberRows=model->numberRows();
    3773       //for (int iRow=0;iRow<numberRows;iRow++)
    3774       //setRowStatus(iRow,ClpSimplex::superBasic);
    3775       std::cout<<"Time before cleanup of full model "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
    3776       primal(1);
    3777       std::cout<<"Total time "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
    3778       delete [] columnCounts;
    3779       delete [] sol;
    3780       delete [] lower;
    3781       delete [] upper;
    3782       delete [] whichBlock;
    3783       delete [] when;
    3784       delete [] columnAdd;
    3785       delete [] rowAdd;
    3786       delete [] elementAdd;
    3787       delete [] objective;
    3788       delete [] top;
    3789       delete [] sub;
    3790     } else {
    3791       delete [] blockInfo;
    3792       printf("What structure - look further?\n");
    3793       loadProblem(*model);
    3794     }
    3795   } else {
    3796     delete [] blockInfo;
    3797     printf("Does not look decomposable????\n");
    3798     loadProblem(*model);
    3799   }
    3800   primal(1);
    3801   return 0;
     3186      if (numberG1<2)
     3187        decomposeType=1;
     3188    }
     3189  }
     3190  if (!decomposeType&&(numberRowBlocks==numberColumnBlocks||
     3191                       numberRowBlocks==numberColumnBlocks-1)) {
     3192    // could be Benders
     3193    int numberG1=0;
     3194    for (int i=0;i<numberColumnBlocks;i++) {
     3195      if (columnCounts[i]>1)
     3196        numberG1++;
     3197    }
     3198    bool masterRows = (numberColumnBlocks==numberRowBlocks);
     3199    if ((masterRows&&numberElementBlocks==2*numberColumnBlocks-1)
     3200        ||(!masterRows&&numberElementBlocks==2*numberColumnBlocks)) {
     3201      if (numberG1<2)
     3202        decomposeType=2;
     3203    }
     3204  }
     3205  delete [] rowCounts;
     3206  delete [] columnCounts;
     3207  delete [] blockInfo;
     3208  // decide what to do
     3209  switch (decomposeType) {
     3210    // No good
     3211  case 0:
     3212    loadProblem(*model,false);
     3213    return dual();
     3214    // DW
     3215  case 1:
     3216    return solveDW(model);
     3217    // Benders
     3218  case 2:
     3219    return solveBenders(model);
     3220  }
     3221  return 0; // to stop compiler warning
    38023222}
    38033223/* This loads a model from a CoinStructuredModel object - returns number of errors.
     
    41033523  return largest;
    41043524}
     3525// Solve using Dantzig-Wolfe decomposition and maybe in parallel
     3526int
     3527ClpSimplex::solveDW(CoinStructuredModel * model)
     3528{
     3529  double time1 = CoinCpuTime();
     3530  int numberColumns = model->numberColumns();
     3531  int numberRowBlocks=model->numberRowBlocks();
     3532  int numberColumnBlocks = model->numberColumnBlocks();
     3533  int numberElementBlocks = model->numberElementBlocks();
     3534  // We already have top level structure
     3535  CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
     3536  for (int i=0;i<numberElementBlocks;i++) {
     3537    CoinModel * thisBlock = model->coinBlock(i);
     3538    assert (thisBlock);
     3539    // just fill in info
     3540    CoinModelBlockInfo info = CoinModelBlockInfo();
     3541    int whatsSet = thisBlock->whatIsSet();
     3542    info.matrix = ((whatsSet&1)!=0) ? 1 : 0;
     3543    info.rhs = ((whatsSet&2)!=0) ? 1 : 0;
     3544    info.rowName = ((whatsSet&4)!=0) ? 1 : 0;
     3545    info.integer = ((whatsSet&32)!=0) ? 1 : 0;
     3546    info.bounds = ((whatsSet&8)!=0) ? 1 : 0;
     3547    info.columnName = ((whatsSet&16)!=0) ? 1 : 0;
     3548    // Which block
     3549    int iRowBlock=model->rowBlock(thisBlock->getRowBlock());
     3550    info.rowBlock=iRowBlock;
     3551    int iColumnBlock=model->columnBlock(thisBlock->getColumnBlock());
     3552    info.columnBlock=iColumnBlock;
     3553    blockInfo[i] = info;
     3554  }
     3555  // make up problems
     3556  int numberBlocks=numberRowBlocks-1;
     3557  // Find master rows and columns
     3558  int * rowCounts = new int [numberRowBlocks];
     3559  CoinZeroN(rowCounts,numberRowBlocks);
     3560  int * columnCounts = new int [numberColumnBlocks+1];
     3561  CoinZeroN(columnCounts,numberColumnBlocks);
     3562  int iBlock;
     3563  for (iBlock = 0;iBlock<numberElementBlocks;iBlock++) {
     3564    int iRowBlock = blockInfo[iBlock].rowBlock;
     3565    rowCounts[iRowBlock]++;
     3566    int iColumnBlock =blockInfo[iBlock].columnBlock;
     3567    columnCounts[iColumnBlock]++;
     3568  }
     3569  int * whichBlock = new int [numberElementBlocks];
     3570  int masterRowBlock=-1;
     3571  for (iBlock = 0;iBlock<numberRowBlocks;iBlock++) {
     3572    if (rowCounts[iBlock]>1) {
     3573      if (masterRowBlock==-1) {
     3574        masterRowBlock=iBlock;
     3575      } else {
     3576        // Can't decode
     3577        masterRowBlock=-2;
     3578        break;
     3579      }
     3580    }
     3581  }
     3582  int masterColumnBlock=-1;
     3583  int kBlock=0;
     3584  for (iBlock = 0;iBlock<numberColumnBlocks;iBlock++) {
     3585    int count=columnCounts[iBlock];
     3586    columnCounts[iBlock]=kBlock;
     3587    kBlock += count;
     3588  }
     3589  for (iBlock = 0;iBlock<numberElementBlocks;iBlock++) {
     3590    int iColumnBlock = blockInfo[iBlock].columnBlock;
     3591    whichBlock[columnCounts[iColumnBlock]]=iBlock;
     3592    columnCounts[iColumnBlock]++;
     3593  }
     3594  for (iBlock = numberColumnBlocks-1;iBlock>=0;iBlock--)
     3595    columnCounts[iBlock+1]=columnCounts[iBlock];
     3596  columnCounts[0]=0;
     3597  for (iBlock = 0;iBlock<numberColumnBlocks;iBlock++) {
     3598    int count=columnCounts[iBlock+1]-columnCounts[iBlock];
     3599    if (count==1) {
     3600      int kBlock = whichBlock[columnCounts[iBlock]];
     3601      int iRowBlock = blockInfo[kBlock].rowBlock;
     3602      if (iRowBlock==masterRowBlock) {
     3603        if (masterColumnBlock==-1) {
     3604          masterColumnBlock=iBlock;
     3605        } else {
     3606          // Can't decode
     3607          masterColumnBlock=-2;
     3608          break;
     3609        }
     3610      }
     3611    }
     3612  }
     3613  if (masterRowBlock<0||masterColumnBlock==-2) {
     3614    // What now
     3615    abort();
     3616  }
     3617  delete [] rowCounts;
     3618  // create all data
     3619  const CoinPackedMatrix ** top = new const CoinPackedMatrix * [numberColumnBlocks];
     3620  ClpSimplex * sub = new ClpSimplex [numberBlocks];
     3621  ClpSimplex master;
     3622  // Set offset
     3623  master.setObjectiveOffset(model->objectiveOffset());
     3624  kBlock=0;
     3625  int masterBlock=-1;
     3626  for (iBlock = 0;iBlock<numberColumnBlocks;iBlock++) {
     3627    top[kBlock]=NULL;
     3628    int start=columnCounts[iBlock];
     3629    int end=columnCounts[iBlock+1];
     3630    assert (end-start<=2);
     3631    for (int j=start;j<end;j++) {
     3632      int jBlock = whichBlock[j];
     3633      int iRowBlock = blockInfo[jBlock].rowBlock;
     3634      int iColumnBlock =blockInfo[jBlock].columnBlock;
     3635      assert (iColumnBlock==iBlock);
     3636      if (iColumnBlock!=masterColumnBlock&&iRowBlock==masterRowBlock) {
     3637        // top matrix
     3638        top[kBlock]=model->coinBlock(jBlock)->packedMatrix();
     3639      } else {
     3640        const CoinPackedMatrix * matrix
     3641          = model->coinBlock(jBlock)->packedMatrix();
     3642        // Get pointers to arrays
     3643        const double * rowLower;
     3644        const double * rowUpper;
     3645        const double * columnLower;
     3646        const double * columnUpper;
     3647        const double * objective;
     3648        model->block(iRowBlock,iColumnBlock,rowLower,rowUpper,
     3649                     columnLower,columnUpper,objective);
     3650        if (iColumnBlock!=masterColumnBlock) {
     3651          // diagonal block
     3652          sub[kBlock].loadProblem(*matrix,columnLower,columnUpper,
     3653                                  objective,rowLower,rowUpper);
     3654          if (true) {
     3655            double * lower = sub[kBlock].columnLower();
     3656            double * upper = sub[kBlock].columnUpper();
     3657            int n = sub[kBlock].numberColumns();
     3658            for (int i=0;i<n;i++) {
     3659              lower[i]=CoinMax(-1.0e8,lower[i]);
     3660              upper[i]=CoinMin(1.0e8,upper[i]);
     3661            }
     3662          }
     3663          if (optimizationDirection_<0.0) {
     3664            double * obj = sub[kBlock].objective();
     3665            int n = sub[kBlock].numberColumns();
     3666            for (int i=0;i<n;i++)
     3667              obj[i] = - obj[i];
     3668          }
     3669          if (this->factorizationFrequency()==200) {
     3670            // User did not touch preset
     3671            sub[kBlock].defaultFactorizationFrequency();
     3672          } else {
     3673            // make sure model has correct value
     3674            sub[kBlock].setFactorizationFrequency(this->factorizationFrequency());
     3675          }
     3676          sub[kBlock].setPerturbation(50);
     3677          // Set columnCounts to be diagonal block index for cleanup
     3678          columnCounts[kBlock]=jBlock;
     3679        } else {
     3680          // master
     3681          masterBlock=jBlock;
     3682          master.loadProblem(*matrix,columnLower,columnUpper,
     3683                             objective,rowLower,rowUpper);
     3684          if (optimizationDirection_<0.0) {
     3685            double * obj = master.objective();
     3686            int n = master.numberColumns();
     3687            for (int i=0;i<n;i++)
     3688              obj[i] = - obj[i];
     3689          }
     3690        }
     3691      }
     3692    }
     3693    if (iBlock!=masterColumnBlock)
     3694      kBlock++;
     3695  }
     3696  delete [] whichBlock;
     3697  delete [] blockInfo;
     3698  // For now master must have been defined (does not have to have columns)
     3699  assert (master.numberRows());
     3700  assert (masterBlock>=0);
     3701  int numberMasterRows = master.numberRows();
     3702  // Overkill in terms of space
     3703  int spaceNeeded = CoinMax(numberBlocks*(numberMasterRows+1),
     3704                            2*numberMasterRows);
     3705  int * rowAdd = new int[spaceNeeded];
     3706  double * elementAdd = new double[spaceNeeded];
     3707  spaceNeeded = numberBlocks;
     3708  int * columnAdd = new int[spaceNeeded+1];
     3709  double * objective = new double[spaceNeeded];
     3710  // Add in costed slacks
     3711  int firstArtificial = master.numberColumns();
     3712  int lastArtificial = firstArtificial;
     3713  if (true) {
     3714    const double * lower = master.rowLower();
     3715    const double * upper = master.rowUpper();
     3716    int kCol=0;
     3717    for (int iRow=0;iRow<numberMasterRows;iRow++) {
     3718      if (lower[iRow]>-1.0e10) {
     3719        rowAdd[kCol]=iRow;
     3720        elementAdd[kCol++]=1.0;
     3721      }
     3722      if (upper[iRow]<1.0e10) {
     3723        rowAdd[kCol]=iRow;
     3724        elementAdd[kCol++]=-1.0;
     3725      }
     3726    }
     3727    if (kCol>spaceNeeded) {
     3728      spaceNeeded = kCol;
     3729      delete [] columnAdd;
     3730      delete [] objective;
     3731      columnAdd = new int[spaceNeeded+1];
     3732      objective = new double[spaceNeeded];
     3733    }
     3734    for (int i=0;i<kCol;i++) {
     3735      columnAdd[i]=i;
     3736      objective[i]=1.0e13;
     3737    }
     3738    columnAdd[kCol]=kCol;
     3739    master.addColumns(kCol,NULL,NULL,objective,
     3740                      columnAdd,rowAdd,elementAdd);
     3741    lastArtificial = master.numberColumns();
     3742  }
     3743  int maxPass=500;
     3744  int iPass;
     3745  double lastObjective=1.0e31;
     3746  // Create convexity rows for proposals
     3747  int numberMasterColumns = master.numberColumns();
     3748  master.resize(numberMasterRows+numberBlocks,numberMasterColumns);
     3749  if (this->factorizationFrequency()==200) {
     3750    // User did not touch preset
     3751    master.defaultFactorizationFrequency();
     3752  } else {
     3753    // make sure model has correct value
     3754    master.setFactorizationFrequency(this->factorizationFrequency());
     3755  }
     3756  master.setPerturbation(50);
     3757  // Arrays to say which block and when created
     3758  int maximumColumns = 2*numberMasterRows+10*numberBlocks;
     3759  whichBlock = new int[maximumColumns];
     3760  int * when = new int[maximumColumns];
     3761  int numberColumnsGenerated=numberBlocks;
     3762  // fill in rhs and add in artificials
     3763  {
     3764    double * rowLower = master.rowLower();
     3765    double * rowUpper = master.rowUpper();
     3766    int iBlock;
     3767    columnAdd[0] = 0;
     3768    for (iBlock=0;iBlock<numberBlocks;iBlock++) {
     3769      int iRow = iBlock + numberMasterRows;;
     3770      rowLower[iRow]=1.0;
     3771      rowUpper[iRow]=1.0;
     3772      rowAdd[iBlock] = iRow;
     3773      elementAdd[iBlock] = 1.0;
     3774      objective[iBlock] = 1.0e13;
     3775      columnAdd[iBlock+1] = iBlock+1;
     3776      when[iBlock]=-1;
     3777      whichBlock[iBlock] = iBlock;
     3778    }
     3779    master.addColumns(numberBlocks,NULL,NULL,objective,
     3780                      columnAdd, rowAdd, elementAdd);
     3781  }
     3782  // and resize matrix to double check clp will be happy
     3783  //master.matrix()->setDimensions(numberMasterRows+numberBlocks,
     3784  //                     numberMasterColumns+numberBlocks);
     3785  std::cout<<"Time to decompose "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
     3786  for (iPass=0;iPass<maxPass;iPass++) {
     3787    printf("Start of pass %d\n",iPass);
     3788    // Solve master - may be infeasible
     3789    //master.scaling(0);
     3790    if (0) {
     3791      master.writeMps("yy.mps");
     3792    }
     3793    // Correct artificials
     3794    double sumArtificials=0.0;
     3795    if (iPass) {
     3796      double * upper = master.columnUpper();
     3797      double * solution = master.primalColumnSolution();
     3798      double * obj = master.objective();
     3799      sumArtificials=0.0;
     3800      for (int i=firstArtificial;i<lastArtificial;i++) {
     3801        sumArtificials += solution[i];
     3802        //assert (solution[i]>-1.0e-2);
     3803        if (solution[i]<1.0e-6) {
     3804#if 0
     3805          // Could take out
     3806          obj[i]=0.0;
     3807          upper[i]=0.0;
     3808#else
     3809          obj[i]=1.0e7;
     3810          upper[i]=1.0e-1;
     3811#endif
     3812          solution[i]=0.0;
     3813          master.setColumnStatus(i,isFixed);
     3814        } else {
     3815          upper[i]=solution[i]+1.0e-5*(1.0+solution[i]);
     3816        }
     3817      }
     3818      printf("Sum of artificials before solve is %g\n",sumArtificials);
     3819    }
     3820    // scale objective to be reasonable
     3821    double scaleFactor = master.scaleObjective(-1.0e9);
     3822    {
     3823      double * dual = master.dualRowSolution();
     3824      int n=master.numberRows();
     3825      memset(dual,0,n*sizeof(double));
     3826      double * solution = master.primalColumnSolution();
     3827      master.clpMatrix()->times(1.0,solution,dual);
     3828      double sum=0.0;
     3829      double * lower = master.rowLower();
     3830      double * upper = master.rowUpper();
     3831      for (int iRow=0;iRow<n;iRow++) {
     3832        double value = dual[iRow];
     3833        if (value>upper[iRow])
     3834          sum += value-upper[iRow];
     3835        else if (value<lower[iRow])
     3836          sum -= value-lower[iRow];
     3837      }
     3838      printf("suminf %g\n",sum);
     3839      lower = master.columnLower();
     3840      upper = master.columnUpper();
     3841      n=master.numberColumns();
     3842      for (int iColumn=0;iColumn<n;iColumn++) {
     3843        double value = solution[iColumn];
     3844        if (value>upper[iColumn]+1.0e-5)
     3845          sum += value-upper[iColumn];
     3846        else if (value<lower[iColumn]-1.0e-5)
     3847          sum -= value-lower[iColumn];
     3848      }
     3849      printf("suminf %g\n",sum);
     3850    }
     3851    master.primal(1);
     3852    // Correct artificials
     3853    sumArtificials=0.0;
     3854    {
     3855      double * solution = master.primalColumnSolution();
     3856      for (int i=firstArtificial;i<lastArtificial;i++) {
     3857        sumArtificials += solution[i];
     3858      }
     3859      printf("Sum of artificials after solve is %g\n",sumArtificials);
     3860    }
     3861    master.scaleObjective(scaleFactor);
     3862    int problemStatus = master.status(); // do here as can change (delcols)
     3863    if (master.numberIterations()==0&&iPass)
     3864      break; // finished
     3865    if (master.objectiveValue()>lastObjective-1.0e-7&&iPass>555)
     3866      break; // finished
     3867    lastObjective = master.objectiveValue();
     3868    // mark basic ones and delete if necessary
     3869    int iColumn;
     3870    numberColumnsGenerated=master.numberColumns()-numberMasterColumns;
     3871    for (iColumn=0;iColumn<numberColumnsGenerated;iColumn++) {
     3872      if (master.getStatus(iColumn+numberMasterColumns)==ClpSimplex::basic)
     3873        when[iColumn]=iPass;
     3874    }
     3875    if (numberColumnsGenerated+numberBlocks>maximumColumns) {
     3876      // delete
     3877      int numberKeep=0;
     3878      int numberDelete=0;
     3879      int * whichDelete = new int[numberColumnsGenerated];
     3880      for (iColumn=0;iColumn<numberColumnsGenerated;iColumn++) {
     3881        if (when[iColumn]>iPass-7) {
     3882          // keep
     3883          when[numberKeep]=when[iColumn];
     3884          whichBlock[numberKeep++]=whichBlock[iColumn];
     3885        } else {
     3886          // delete
     3887          whichDelete[numberDelete++]=iColumn+numberMasterColumns;
     3888        }
     3889      }
     3890      numberColumnsGenerated -= numberDelete;
     3891      master.deleteColumns(numberDelete,whichDelete);
     3892      delete [] whichDelete;
     3893    }
     3894    const double * dual=NULL;
     3895    bool deleteDual=false;
     3896    if (problemStatus==0) {
     3897      dual = master.dualRowSolution();
     3898    } else if (problemStatus==1) {
     3899      // could do composite objective
     3900      dual = master.infeasibilityRay();
     3901      deleteDual = true;
     3902      printf("The sum of infeasibilities is %g\n",
     3903             master.sumPrimalInfeasibilities());
     3904    } else if (!master.numberColumns()) {
     3905      assert(!iPass);
     3906      dual = master.dualRowSolution();
     3907      memset(master.dualRowSolution(),
     3908             0,(numberMasterRows+numberBlocks)*sizeof(double));
     3909    } else {
     3910      abort();
     3911    }
     3912    // Scale back on first time
     3913    if (!iPass) {
     3914      double * dual2 = master.dualRowSolution();
     3915      for (int i=0;i<numberMasterRows+numberBlocks;i++) {
     3916        dual2[i] *= 1.0e-7;
     3917      }
     3918      dual = master.dualRowSolution();
     3919    }
     3920    // Create objective for sub problems and solve
     3921    columnAdd[0]=0;
     3922    int numberProposals=0;
     3923    for (iBlock=0;iBlock<numberBlocks;iBlock++) {
     3924      int numberColumns2 = sub[iBlock].numberColumns();
     3925      double * saveObj = new double [numberColumns2];
     3926      double * objective2 = sub[iBlock].objective();
     3927      memcpy(saveObj,objective2,numberColumns2*sizeof(double));
     3928      // new objective
     3929      top[iBlock]->transposeTimes(dual,objective2);
     3930      int i;
     3931      if (problemStatus==0) {
     3932        for (i=0;i<numberColumns2;i++)
     3933          objective2[i] = saveObj[i]-objective2[i];
     3934      } else {
     3935        for (i=0;i<numberColumns2;i++)
     3936          objective2[i] = -objective2[i];
     3937      }
     3938      // scale objective to be reasonable
     3939      double scaleFactor =
     3940        sub[iBlock].scaleObjective((sumArtificials>1.0e-5) ? -1.0e-4 :-1.0e9);
     3941      if (iPass) {
     3942        sub[iBlock].primal();
     3943      } else {
     3944        sub[iBlock].dual();
     3945      }
     3946      sub[iBlock].scaleObjective(scaleFactor);
     3947      if (!sub[iBlock].isProvenOptimal()&&
     3948          !sub[iBlock].isProvenDualInfeasible()) {
     3949        memset(objective2,0,numberColumns2*sizeof(double));
     3950        sub[iBlock].primal();
     3951        if (problemStatus==0) {
     3952          for (i=0;i<numberColumns2;i++)
     3953            objective2[i] = saveObj[i]-objective2[i];
     3954        } else {
     3955          for (i=0;i<numberColumns2;i++)
     3956            objective2[i] = -objective2[i];
     3957        }
     3958        double scaleFactor = sub[iBlock].scaleObjective(-1.0e9);
     3959        sub[iBlock].primal(1);
     3960        sub[iBlock].scaleObjective(scaleFactor);
     3961      }
     3962      memcpy(objective2,saveObj,numberColumns2*sizeof(double));
     3963      // get proposal
     3964      if (sub[iBlock].numberIterations()||!iPass) {
     3965        double objValue =0.0;
     3966        int start = columnAdd[numberProposals];
     3967        // proposal
     3968        if (sub[iBlock].isProvenOptimal()) {
     3969          const double * solution = sub[iBlock].primalColumnSolution();
     3970          top[iBlock]->times(solution,elementAdd+start);
     3971          for (i=0;i<numberColumns2;i++)
     3972            objValue += solution[i]*saveObj[i];
     3973          // See if good dj and pack down
     3974          int number=start;
     3975          double dj = objValue;
     3976          if (problemStatus)
     3977            dj=0.0;
     3978          double smallest=1.0e100;
     3979          double largest=0.0;
     3980          for (i=0;i<numberMasterRows;i++) {
     3981            double value = elementAdd[start+i];
     3982            if (fabs(value)>1.0e-15) {
     3983              dj -= dual[i]*value;
     3984              smallest = CoinMin(smallest,fabs(value));
     3985              largest = CoinMax(largest,fabs(value));
     3986              rowAdd[number]=i;
     3987              elementAdd[number++]=value;
     3988            }
     3989          }
     3990          // and convexity
     3991          dj -= dual[numberMasterRows+iBlock];
     3992          rowAdd[number]=numberMasterRows+iBlock;
     3993          elementAdd[number++]=1.0;
     3994          // if elements large then scale?
     3995          //if (largest>1.0e8||smallest<1.0e-8)
     3996          printf("For subproblem %d smallest - %g, largest %g - dj %g\n",
     3997                 iBlock,smallest,largest,dj);
     3998          if (dj<-1.0e-6||!iPass) {
     3999            // take
     4000            objective[numberProposals]=objValue;
     4001            columnAdd[++numberProposals]=number;
     4002            when[numberColumnsGenerated]=iPass;
     4003            whichBlock[numberColumnsGenerated++]=iBlock;
     4004          }
     4005        } else if (sub[iBlock].isProvenDualInfeasible()) {
     4006          // use ray
     4007          const double * solution = sub[iBlock].unboundedRay();
     4008          top[iBlock]->times(solution,elementAdd+start);
     4009          for (i=0;i<numberColumns2;i++)
     4010            objValue += solution[i]*saveObj[i];
     4011          // See if good dj and pack down
     4012          int number=start;
     4013          double dj = objValue;
     4014          double smallest=1.0e100;
     4015          double largest=0.0;
     4016          for (i=0;i<numberMasterRows;i++) {
     4017            double value = elementAdd[start+i];
     4018            if (fabs(value)>1.0e-15) {
     4019              dj -= dual[i]*value;
     4020              smallest = CoinMin(smallest,fabs(value));
     4021              largest = CoinMax(largest,fabs(value));
     4022              rowAdd[number]=i;
     4023              elementAdd[number++]=value;
     4024            }
     4025          }
     4026          // if elements large or small then scale?
     4027          //if (largest>1.0e8||smallest<1.0e-8)
     4028          printf("For subproblem ray %d smallest - %g, largest %g - dj %g\n",
     4029                 iBlock,smallest,largest,dj);
     4030          if (dj<-1.0e-6) {
     4031            // take
     4032            objective[numberProposals]=objValue;
     4033            columnAdd[++numberProposals]=number;
     4034            when[numberColumnsGenerated]=iPass;
     4035            whichBlock[numberColumnsGenerated++]=iBlock;
     4036          }
     4037        } else {
     4038          abort();
     4039        }
     4040      }
     4041      delete [] saveObj;
     4042    }
     4043    if (deleteDual)
     4044      delete [] dual;
     4045    if (numberProposals)
     4046      master.addColumns(numberProposals,NULL,NULL,objective,
     4047                        columnAdd,rowAdd,elementAdd);
     4048  }
     4049  std::cout<<"Time at end of D-W "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
     4050  //master.scaling(0);
     4051  //master.primal(1);
     4052  loadProblem(*model);
     4053  // now put back a good solution
     4054  double * lower = new double[numberMasterRows+numberBlocks];
     4055  double * upper = new double[numberMasterRows+numberBlocks];
     4056  numberColumnsGenerated  += numberMasterColumns;
     4057  double * sol = new double[numberColumnsGenerated];
     4058  const double * solution = master.primalColumnSolution();
     4059  const double * masterLower = master.rowLower();
     4060  const double * masterUpper = master.rowUpper();
     4061  double * fullSolution = primalColumnSolution();
     4062  const double * fullLower = columnLower();
     4063  const double * fullUpper = columnUpper();
     4064  const double * rowSolution = master.primalRowSolution();
     4065  double * fullRowSolution = primalRowSolution();
     4066  const int * rowBack = model->coinBlock(masterBlock)->originalRows();
     4067  int numberRows2 = model->coinBlock(masterBlock)->numberRows();
     4068  const int * columnBack = model->coinBlock(masterBlock)->originalColumns();
     4069  int numberColumns2 = model->coinBlock(masterBlock)->numberColumns();
     4070  for (int iRow=0;iRow<numberRows2;iRow++) {
     4071    int kRow = rowBack[iRow];
     4072    setRowStatus(kRow,master.getRowStatus(iRow));
     4073    fullRowSolution[kRow]=rowSolution[iRow];
     4074  }
     4075  for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
     4076    int kColumn = columnBack[iColumn];
     4077    setStatus(kColumn,master.getStatus(iColumn));
     4078    fullSolution[kColumn]=solution[iColumn];
     4079  }
     4080  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
     4081    // move basis
     4082    int kBlock = columnCounts[iBlock];
     4083    const int * rowBack = model->coinBlock(kBlock)->originalRows();
     4084    int numberRows2 = model->coinBlock(kBlock)->numberRows();
     4085    const int * columnBack = model->coinBlock(kBlock)->originalColumns();
     4086    int numberColumns2 = model->coinBlock(kBlock)->numberColumns();
     4087    for (int iRow=0;iRow<numberRows2;iRow++) {
     4088      int kRow = rowBack[iRow];
     4089      setRowStatus(kRow,sub[iBlock].getRowStatus(iRow));
     4090    }
     4091    for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
     4092      int kColumn = columnBack[iColumn];
     4093      setStatus(kColumn,sub[iBlock].getStatus(iColumn));
     4094    }
     4095    // convert top bit to by rows
     4096    CoinPackedMatrix topMatrix = *top[iBlock];
     4097    topMatrix.reverseOrdering();
     4098    // zero solution
     4099    memset(sol,0,numberColumnsGenerated*sizeof(double));
     4100   
     4101    for (int i=numberMasterColumns;i<numberColumnsGenerated;i++) {
     4102      if (whichBlock[i-numberMasterColumns]==iBlock)
     4103        sol[i] = solution[i];
     4104    }
     4105    memset(lower,0,(numberMasterRows+numberBlocks)*sizeof(double));
     4106    master.clpMatrix()->times(1.0,sol,lower);
     4107    for (int iRow=0;iRow<numberMasterRows;iRow++) {
     4108      double value = lower[iRow];
     4109      if (masterUpper[iRow]<1.0e20)
     4110        upper[iRow] = value;
     4111      else
     4112        upper[iRow]=COIN_DBL_MAX;
     4113      if (masterLower[iRow]>-1.0e20)
     4114        lower[iRow] = value;
     4115      else
     4116        lower[iRow]=-COIN_DBL_MAX;
     4117    }
     4118    sub[iBlock].addRows(numberMasterRows,lower,upper,
     4119                        topMatrix.getVectorStarts(),
     4120                        topMatrix.getVectorLengths(),
     4121                        topMatrix.getIndices(),
     4122                        topMatrix.getElements());
     4123    sub[iBlock].primal(1);
     4124    const double * subSolution = sub[iBlock].primalColumnSolution();
     4125    const double * subRowSolution = sub[iBlock].primalRowSolution();
     4126    // move solution
     4127    for (int iRow=0;iRow<numberRows2;iRow++) {
     4128      int kRow = rowBack[iRow];
     4129      fullRowSolution[kRow]=subRowSolution[iRow];
     4130    }
     4131    for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
     4132      int kColumn = columnBack[iColumn];
     4133      fullSolution[kColumn]=subSolution[iColumn];
     4134    }
     4135  }
     4136  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     4137    if (fullSolution[iColumn]<fullUpper[iColumn]-1.0e-8&&
     4138        fullSolution[iColumn]>fullLower[iColumn]+1.0e-8) {
     4139      if (getStatus(iColumn)!=ClpSimplex::basic) {
     4140        if (columnLower_[iColumn]>-1.0e30||
     4141            columnUpper_[iColumn]<1.0e30)
     4142          setStatus(iColumn,ClpSimplex::superBasic);
     4143        else
     4144          setStatus(iColumn,ClpSimplex::isFree);
     4145      }
     4146    } else if (fullSolution[iColumn]>=fullUpper[iColumn]-1.0e-8) {
     4147      // may help to make rest non basic
     4148      if (getStatus(iColumn)!=ClpSimplex::basic)
     4149        setStatus(iColumn,ClpSimplex::atUpperBound);
     4150    } else if (fullSolution[iColumn]<=fullLower[iColumn]+1.0e-8) {
     4151      // may help to make rest non basic
     4152      if (getStatus(iColumn)!=ClpSimplex::basic)
     4153        setStatus(iColumn,ClpSimplex::atLowerBound);
     4154    }
     4155  }
     4156  //int numberRows=model->numberRows();
     4157  //for (int iRow=0;iRow<numberRows;iRow++)
     4158  //setRowStatus(iRow,ClpSimplex::superBasic);
     4159  std::cout<<"Time before cleanup of full model "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
     4160  primal(1);
     4161  std::cout<<"Total time "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
     4162  delete [] columnCounts;
     4163  delete [] sol;
     4164  delete [] lower;
     4165  delete [] upper;
     4166  delete [] whichBlock;
     4167  delete [] when;
     4168  delete [] columnAdd;
     4169  delete [] rowAdd;
     4170  delete [] elementAdd;
     4171  delete [] objective;
     4172  delete [] top;
     4173  delete [] sub;
     4174  return 0;
     4175}
     4176// Solve using Benders decomposition and maybe in parallel
     4177int
     4178ClpSimplex::solveBenders(CoinStructuredModel * model)
     4179{
     4180  double time1 = CoinCpuTime();
     4181  //int numberColumns = model->numberColumns();
     4182  int numberRowBlocks=model->numberRowBlocks();
     4183  int numberColumnBlocks = model->numberColumnBlocks();
     4184  int numberElementBlocks = model->numberElementBlocks();
     4185  // We already have top level structure
     4186  CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
     4187  for (int i=0;i<numberElementBlocks;i++) {
     4188    CoinModel * thisBlock = model->coinBlock(i);
     4189    assert (thisBlock);
     4190    // just fill in info
     4191    CoinModelBlockInfo info = CoinModelBlockInfo();
     4192    int whatsSet = thisBlock->whatIsSet();
     4193    info.matrix = ((whatsSet&1)!=0) ? 1 : 0;
     4194    info.rhs = ((whatsSet&2)!=0) ? 1 : 0;
     4195    info.rowName = ((whatsSet&4)!=0) ? 1 : 0;
     4196    info.integer = ((whatsSet&32)!=0) ? 1 : 0;
     4197    info.bounds = ((whatsSet&8)!=0) ? 1 : 0;
     4198    info.columnName = ((whatsSet&16)!=0) ? 1 : 0;
     4199    // Which block
     4200    int iRowBlock=model->rowBlock(thisBlock->getRowBlock());
     4201    info.rowBlock=iRowBlock;
     4202    int iColumnBlock=model->columnBlock(thisBlock->getColumnBlock());
     4203    info.columnBlock=iColumnBlock;
     4204    blockInfo[i] = info;
     4205  }
     4206  // make up problems
     4207  int numberBlocks=numberColumnBlocks-1;
     4208  // Find master columns and rows
     4209  int * columnCounts = new int [numberColumnBlocks];
     4210  CoinZeroN(columnCounts,numberColumnBlocks);
     4211  int * rowCounts = new int [numberRowBlocks+1];
     4212  CoinZeroN(rowCounts,numberRowBlocks);
     4213  int iBlock;
     4214  for (iBlock = 0;iBlock<numberElementBlocks;iBlock++) {
     4215    int iColumnBlock = blockInfo[iBlock].columnBlock;
     4216    columnCounts[iColumnBlock]++;
     4217    int iRowBlock =blockInfo[iBlock].rowBlock;
     4218    rowCounts[iRowBlock]++;
     4219  }
     4220  int * whichBlock = new int [numberElementBlocks];
     4221  int masterColumnBlock=-1;
     4222  for (iBlock = 0;iBlock<numberColumnBlocks;iBlock++) {
     4223    if (columnCounts[iBlock]>1) {
     4224      if (masterColumnBlock==-1) {
     4225        masterColumnBlock=iBlock;
     4226      } else {
     4227        // Can't decode
     4228        masterColumnBlock=-2;
     4229        break;
     4230      }
     4231    }
     4232  }
     4233  int masterRowBlock=-1;
     4234  int kBlock=0;
     4235  for (iBlock = 0;iBlock<numberRowBlocks;iBlock++) {
     4236    int count=rowCounts[iBlock];
     4237    rowCounts[iBlock]=kBlock;
     4238    kBlock += count;
     4239  }
     4240  for (iBlock = 0;iBlock<numberElementBlocks;iBlock++) {
     4241    int iRowBlock = blockInfo[iBlock].rowBlock;
     4242    whichBlock[rowCounts[iRowBlock]]=iBlock;
     4243    rowCounts[iRowBlock]++;
     4244  }
     4245  for (iBlock = numberRowBlocks-1;iBlock>=0;iBlock--)
     4246    rowCounts[iBlock+1]=rowCounts[iBlock];
     4247  rowCounts[0]=0;
     4248  for (iBlock = 0;iBlock<numberRowBlocks;iBlock++) {
     4249    int count=rowCounts[iBlock+1]-rowCounts[iBlock];
     4250    if (count==1) {
     4251      int kBlock = whichBlock[rowCounts[iBlock]];
     4252      int iColumnBlock = blockInfo[kBlock].columnBlock;
     4253      if (iColumnBlock==masterColumnBlock) {
     4254        if (masterRowBlock==-1) {
     4255          masterRowBlock=iBlock;
     4256        } else {
     4257          // Can't decode
     4258          masterRowBlock=-2;
     4259          break;
     4260        }
     4261      }
     4262    }
     4263  }
     4264  if (masterColumnBlock<0||masterRowBlock==-2) {
     4265    // What now
     4266    abort();
     4267  }
     4268  delete [] columnCounts;
     4269  // create all data
     4270  const CoinPackedMatrix ** first = new const CoinPackedMatrix * [numberRowBlocks];
     4271  ClpSimplex * sub = new ClpSimplex [numberBlocks];
     4272  ClpSimplex master;
     4273  // Set offset
     4274  master.setObjectiveOffset(model->objectiveOffset());
     4275  kBlock=0;
     4276  int masterBlock=-1;
     4277  for (iBlock = 0;iBlock<numberRowBlocks;iBlock++) {
     4278    first[kBlock]=NULL;
     4279    int start=rowCounts[iBlock];
     4280    int end=rowCounts[iBlock+1];
     4281    assert (end-start<=2);
     4282    for (int j=start;j<end;j++) {
     4283      int jBlock = whichBlock[j];
     4284      int iColumnBlock = blockInfo[jBlock].columnBlock;
     4285      int iRowBlock =blockInfo[jBlock].rowBlock;
     4286      assert (iRowBlock==iBlock);
     4287      if (iRowBlock!=masterRowBlock&&iColumnBlock==masterColumnBlock) {
     4288        // first matrix
     4289        first[kBlock]=model->coinBlock(jBlock)->packedMatrix();
     4290      } else {
     4291        const CoinPackedMatrix * matrix
     4292          = model->coinBlock(jBlock)->packedMatrix();
     4293        // Get pointers to arrays
     4294        const double * columnLower;
     4295        const double * columnUpper;
     4296        const double * rowLower;
     4297        const double * rowUpper;
     4298        const double * objective;
     4299        model->block(iRowBlock,iColumnBlock,rowLower,rowUpper,
     4300                     columnLower,columnUpper,objective);
     4301        if (iRowBlock!=masterRowBlock) {
     4302          // diagonal block
     4303          sub[kBlock].loadProblem(*matrix,columnLower,columnUpper,
     4304                                  objective,rowLower,rowUpper);
     4305          if (optimizationDirection_<0.0) {
     4306            double * obj = sub[kBlock].objective();
     4307            int n = sub[kBlock].numberColumns();
     4308            for (int i=0;i<n;i++)
     4309              obj[i] = - obj[i];
     4310          }
     4311          if (this->factorizationFrequency()==200) {
     4312            // User did not touch preset
     4313            sub[kBlock].defaultFactorizationFrequency();
     4314          } else {
     4315            // make sure model has correct value
     4316            sub[kBlock].setFactorizationFrequency(this->factorizationFrequency());
     4317          }
     4318          sub[kBlock].setPerturbation(50);
     4319          // Set rowCounts to be diagonal block index for cleanup
     4320          rowCounts[kBlock]=jBlock;
     4321        } else {
     4322          // master
     4323          masterBlock=jBlock;
     4324          master.loadProblem(*matrix,columnLower,columnUpper,
     4325                             objective,rowLower,rowUpper);
     4326          if (optimizationDirection_<0.0) {
     4327            double * obj = master.objective();
     4328            int n = master.numberColumns();
     4329            for (int i=0;i<n;i++)
     4330              obj[i] = - obj[i];
     4331          }
     4332        }
     4333      }
     4334    }
     4335    if (iBlock!=masterRowBlock)
     4336      kBlock++;
     4337  }
     4338  delete [] whichBlock;
     4339  delete [] blockInfo;
     4340  // For now master must have been defined (does not have to have rows)
     4341  assert (master.numberColumns());
     4342  assert (masterBlock>=0);
     4343  int numberMasterColumns = master.numberColumns();
     4344  // Overkill in terms of space
     4345  int spaceNeeded = CoinMax(numberBlocks*(numberMasterColumns+1),
     4346                            2*numberMasterColumns);
     4347  int * columnAdd = new int[spaceNeeded];
     4348  double * elementAdd = new double[spaceNeeded];
     4349  spaceNeeded = numberBlocks;
     4350  int * rowAdd = new int[spaceNeeded+1];
     4351  double * objective = new double[spaceNeeded];
     4352  int maxPass=500;
     4353  int iPass;
     4354  double lastObjective=-1.0e31;
     4355  // Create columns for proposals
     4356  int numberMasterRows = master.numberRows();
     4357  master.resize(numberMasterColumns+numberBlocks,numberMasterRows);
     4358  if (this->factorizationFrequency()==200) {
     4359    // User did not touch preset
     4360    master.defaultFactorizationFrequency();
     4361  } else {
     4362    // make sure model has correct value
     4363    master.setFactorizationFrequency(this->factorizationFrequency());
     4364  }
     4365  master.setPerturbation(50);
     4366  // Arrays to say which block and when created
     4367  int maximumRows = 2*numberMasterColumns+10*numberBlocks;
     4368  whichBlock = new int[maximumRows];
     4369  int * when = new int[maximumRows];
     4370  int numberRowsGenerated=numberBlocks;
     4371  // Add extra variables
     4372  {
     4373    int iBlock;
     4374    columnAdd[0] = 0;
     4375    for (iBlock=0;iBlock<numberBlocks;iBlock++) {
     4376      objective[iBlock] = 1.0;
     4377      columnAdd[iBlock+1] = 0;
     4378      when[iBlock]=-1;
     4379      whichBlock[iBlock] = iBlock;
     4380    }
     4381    master.addColumns(numberBlocks,NULL,NULL,objective,
     4382                      columnAdd, rowAdd, elementAdd);
     4383  }
     4384  std::cout<<"Time to decompose "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
     4385  for (iPass=0;iPass<maxPass;iPass++) {
     4386    printf("Start of pass %d\n",iPass);
     4387    // Solve master - may be unbounded
     4388    //master.scaling(0);
     4389    if (1) {
     4390      master.writeMps("yy.mps");
     4391    }
     4392    master.dual();
     4393    int problemStatus = master.status(); // do here as can change (delcols)
     4394    if (master.numberIterations()==0&&iPass)
     4395      break; // finished
     4396    if (master.objectiveValue()<lastObjective+1.0e-7&&iPass>555)
     4397      break; // finished
     4398    lastObjective = master.objectiveValue();
     4399    // mark non-basic rows and delete if necessary
     4400    int iRow;
     4401    numberRowsGenerated=master.numberRows()-numberMasterRows;
     4402    for (iRow=0;iRow<numberRowsGenerated;iRow++) {
     4403      if (master.getStatus(iRow+numberMasterRows)!=ClpSimplex::basic)
     4404        when[iRow]=iPass;
     4405    }
     4406    if (numberRowsGenerated>maximumRows) {
     4407      // delete
     4408      int numberKeep=0;
     4409      int numberDelete=0;
     4410      int * whichDelete = new int[numberRowsGenerated];
     4411      for (iRow=0;iRow<numberRowsGenerated;iRow++) {
     4412        if (when[iRow]>iPass-7) {
     4413          // keep
     4414          when[numberKeep]=when[iRow];
     4415          whichBlock[numberKeep++]=whichBlock[iRow];
     4416        } else {
     4417          // delete
     4418          whichDelete[numberDelete++]=iRow+numberMasterRows;
     4419        }
     4420      }
     4421      numberRowsGenerated -= numberDelete;
     4422      master.deleteRows(numberDelete,whichDelete);
     4423      delete [] whichDelete;
     4424    }
     4425    const double * primal=NULL;
     4426    bool deletePrimal=false;
     4427    if (problemStatus==0) {
     4428      primal = master.primalColumnSolution();
     4429    } else if (problemStatus==2) {
     4430      // could do composite objective
     4431      primal = master.unboundedRay();
     4432      deletePrimal = true;
     4433      printf("The sum of infeasibilities is %g\n",
     4434             master.sumPrimalInfeasibilities());
     4435    } else if (!master.numberRows()) {
     4436      assert(!iPass);
     4437      primal = master.primalColumnSolution();
     4438      memset(master.primalColumnSolution(),
     4439             0,numberMasterColumns*sizeof(double));
     4440    } else {
     4441      abort();
     4442    }
     4443    // Create rhs for sub problems and solve
     4444    rowAdd[0]=0;
     4445    int numberProposals=0;
     4446    for (iBlock=0;iBlock<numberBlocks;iBlock++) {
     4447      int numberRows2 = sub[iBlock].numberRows();
     4448      double * saveLower = new double [numberRows2];
     4449      double * lower2 = sub[iBlock].rowLower();
     4450      double * saveUpper = new double [numberRows2];
     4451      double * upper2 = sub[iBlock].rowUpper();
     4452      // new rhs
     4453      CoinZeroN(saveUpper,numberRows2);
     4454      first[iBlock]->times(primal,saveUpper);
     4455      for (int i=0;i<numberRows2;i++) {
     4456        double value = saveUpper[i];
     4457        saveLower[i]=lower2[i];
     4458        saveUpper[i]=upper2[i];
     4459        if (saveLower[i]>-1.0e30)
     4460          lower2[i] -= value;
     4461        if (saveUpper[i]<1.0e30)
     4462          upper2[i] -= value;
     4463      }
     4464      sub[iBlock].dual();
     4465      memcpy(lower2,saveLower,numberRows2*sizeof(double));
     4466      memcpy(upper2,saveUpper,numberRows2*sizeof(double));
     4467      // get proposal
     4468      if (sub[iBlock].numberIterations()||!iPass) {
     4469        double objValue =0.0;
     4470        int start = rowAdd[numberProposals];
     4471        // proposal
     4472        if (sub[iBlock].isProvenOptimal()) {
     4473          const double * solution = sub[iBlock].dualRowSolution();
     4474          first[iBlock]->transposeTimes(solution,elementAdd+start);
     4475          for (int i=0;i<numberRows2;i++) {
     4476            if (solution[i]<-dualTolerance_) {
     4477              // at upper
     4478              assert (saveUpper[i]<1.0e30);
     4479              objValue += solution[i]*saveUpper[i];
     4480            } else if (solution[i]>dualTolerance_) {
     4481              // at lower
     4482              assert (saveLower[i]>-1.0e30);
     4483              objValue += solution[i]*saveLower[i];
     4484            }
     4485          }
     4486         
     4487          // See if cuts off and pack down
     4488          int number=start;
     4489          double infeas = objValue;
     4490          double smallest=1.0e100;
     4491          double largest=0.0;
     4492          for (int i=0;i<numberMasterColumns;i++) {
     4493            double value = elementAdd[start+i];
     4494            if (fabs(value)>1.0e-15) {
     4495              infeas -= primal[i]*value;
     4496              smallest = CoinMin(smallest,fabs(value));
     4497              largest = CoinMax(largest,fabs(value));
     4498              columnAdd[number]=i;
     4499              elementAdd[number++]=-value;
     4500            }
     4501          }
     4502          columnAdd[number]=numberMasterColumns+iBlock;
     4503          elementAdd[number++]=-1.0;
     4504          // if elements large then scale?
     4505          //if (largest>1.0e8||smallest<1.0e-8)
     4506          printf("For subproblem %d smallest - %g, largest %g - infeas %g\n",
     4507                 iBlock,smallest,largest,infeas);
     4508          if (infeas<-1.0e-6||!iPass) {
     4509            // take
     4510            objective[numberProposals]=objValue;
     4511            rowAdd[++numberProposals]=number;
     4512            when[numberRowsGenerated]=iPass;
     4513            whichBlock[numberRowsGenerated++]=iBlock;
     4514          }
     4515        } else if (sub[iBlock].isProvenPrimalInfeasible()) {
     4516          // use ray
     4517          const double * solution = sub[iBlock].infeasibilityRay();
     4518          first[iBlock]->transposeTimes(solution,elementAdd+start);
     4519          for (int i=0;i<numberRows2;i++) {
     4520            if (solution[i]<-dualTolerance_) {
     4521              // at upper
     4522              assert (saveUpper[i]<1.0e30);
     4523              objValue += solution[i]*saveUpper[i];
     4524            } else if (solution[i]>dualTolerance_) {
     4525              // at lower
     4526              assert (saveLower[i]>-1.0e30);
     4527              objValue += solution[i]*saveLower[i];
     4528            }
     4529          }
     4530          // See if good infeas and pack down
     4531          int number=start;
     4532          double infeas = objValue;
     4533          double smallest=1.0e100;
     4534          double largest=0.0;
     4535          for (int i=0;i<numberMasterColumns;i++) {
     4536            double value = elementAdd[start+i];
     4537            if (fabs(value)>1.0e-15) {
     4538              infeas -= primal[i]*value;
     4539              smallest = CoinMin(smallest,fabs(value));
     4540              largest = CoinMax(largest,fabs(value));
     4541              columnAdd[number]=i;
     4542              elementAdd[number++]=-value;
     4543            }
     4544          }
     4545          // if elements large or small then scale?
     4546          //if (largest>1.0e8||smallest<1.0e-8)
     4547          printf("For subproblem ray %d smallest - %g, largest %g - infeas %g\n",
     4548                 iBlock,smallest,largest,infeas);
     4549          if (infeas<-1.0e-6) {
     4550            // take
     4551            objective[numberProposals]=objValue;
     4552            rowAdd[++numberProposals]=number;
     4553            when[numberRowsGenerated]=iPass;
     4554            whichBlock[numberRowsGenerated++]=iBlock;
     4555          }
     4556        } else {
     4557          abort();
     4558        }
     4559      }
     4560      delete [] saveLower;
     4561      delete [] saveUpper;
     4562    }
     4563    if (deletePrimal)
     4564      delete [] primal;
     4565    if (numberProposals) {
     4566      master.addRows(numberProposals,NULL,objective,
     4567                        rowAdd,columnAdd,elementAdd);
     4568    }
     4569  }
     4570  std::cout<<"Time at end of Benders "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
     4571  //master.scaling(0);
     4572  //master.primal(1);
     4573  loadProblem(*model);
     4574  // now put back a good solution
     4575  const double * columnSolution = master.primalColumnSolution();
     4576  double * fullColumnSolution = primalColumnSolution();
     4577  const int * columnBack = model->coinBlock(masterBlock)->originalColumns();
     4578  int numberColumns2 = model->coinBlock(masterBlock)->numberColumns();
     4579  const int * rowBack = model->coinBlock(masterBlock)->originalRows();
     4580  int numberRows2 = model->coinBlock(masterBlock)->numberRows();
     4581  for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
     4582    int kColumn = columnBack[iColumn];
     4583    setColumnStatus(kColumn,master.getColumnStatus(iColumn));
     4584    fullColumnSolution[kColumn]=columnSolution[iColumn];
     4585  }
     4586  for (int iRow=0;iRow<numberRows2;iRow++) {
     4587    int kRow = rowBack[iRow];
     4588    setStatus(kRow,master.getStatus(iRow));
     4589    //fullSolution[kRow]=solution[iRow];
     4590  }
     4591  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
     4592    // move basis
     4593    int kBlock = rowCounts[iBlock];
     4594    const int * columnBack = model->coinBlock(kBlock)->originalColumns();
     4595    int numberColumns2 = model->coinBlock(kBlock)->numberColumns();
     4596    const int * rowBack = model->coinBlock(kBlock)->originalRows();
     4597    int numberRows2 = model->coinBlock(kBlock)->numberRows();
     4598    const double * subColumnSolution = sub[iBlock].primalColumnSolution();
     4599    for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
     4600      int kColumn = columnBack[iColumn];
     4601      setColumnStatus(kColumn,sub[iBlock].getColumnStatus(iColumn));
     4602      fullColumnSolution[kColumn]=subColumnSolution[iColumn];
     4603    }
     4604    for (int iRow=0;iRow<numberRows2;iRow++) {
     4605      int kRow = rowBack[iRow];
     4606      setStatus(kRow,sub[iBlock].getStatus(iRow));
     4607      setStatus(kRow,atLowerBound);
     4608    }
     4609  }
     4610  double * fullSolution = primalRowSolution();
     4611  CoinZeroN(fullSolution,numberRows_);
     4612  times(1.0,fullColumnSolution,fullSolution);
     4613  //int numberColumns=model->numberColumns();
     4614  //for (int iColumn=0;iColumn<numberColumns;iColumn++)
     4615  //setColumnStatus(iColumn,ClpSimplex::superBasic);
     4616  std::cout<<"Time before cleanup of full model "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
     4617  this->primal(1);
     4618  std::cout<<"Total time "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
     4619  delete [] rowCounts;
     4620  //delete [] sol;
     4621  //delete [] lower;
     4622  //delete [] upper;
     4623  delete [] whichBlock;
     4624  delete [] when;
     4625  delete [] rowAdd;
     4626  delete [] columnAdd;
     4627  delete [] elementAdd;
     4628  delete [] objective;
     4629  delete [] first;
     4630  delete [] sub;
     4631  return 0;
     4632}
Note: See TracChangeset for help on using the changeset viewer.