Changeset 1326
- Timestamp:
- Jan 20, 2009 2:19:35 PM (11 years ago)
- Location:
- trunk/Clp
- Files:
-
- 1 added
- 3 edited
-
configure (modified) (15 diffs)
-
examples/decomp3.cpp (added)
-
src/ClpSimplex.hpp (modified) (1 diff)
-
src/ClpSolve.cpp (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/Clp/configure
r1264 r1326 3505 3505 ac_compiler_gnu=$ac_cv_cxx_compiler_gnu 3506 3506 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 3510 ac_ext=cc 3511 ac_cpp='$CXXCPP $CPPFLAGS' 3512 ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5' 3513 ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' 3514 ac_compiler_gnu=$ac_cv_cxx_compiler_gnu 3515 3516 echo "$as_me:$LINENO: checking whether C++ compiler $CXX works" >&5 3517 echo $ECHO_N "checking whether C++ compiler $CXX works... $ECHO_C" >&6; 3518 cat >conftest.$ac_ext <<_ACEOF 3519 /* confdefs.h. */ 3520 _ACEOF 3521 cat confdefs.h >>conftest.$ac_ext 3522 cat >>conftest.$ac_ext <<_ACEOF 3523 /* end confdefs.h. */ 3524 3525 int 3526 main () 3527 { 3528 int i=0; 3529 ; 3530 return 0; 3531 } 3532 _ACEOF 3533 rm -f conftest.$ac_objext 3534 if { (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 3556 echo "${ECHO_T}yes" >&6 3557 else 3558 echo "$as_me: failed program was:" >&5 3559 sed 's/^/| /' conftest.$ac_ext >&5 3560 3561 echo "$as_me:$LINENO: result: no" >&5 3562 echo "${ECHO_T}no" >&6 3563 { { echo "$as_me:$LINENO: error: failed to find a C++ compiler or C++ compiler $CXX does not work" >&5 3564 echo "$as_me: error: failed to find a C++ compiler or C++ compiler $CXX does not work" >&2;} 3510 3565 { (exit 1); exit 1; }; } 3511 fi 3566 3567 fi 3568 rm -f conftest.err conftest.$ac_objext conftest.$ac_ext 3569 ac_ext=cc 3570 ac_cpp='$CXXCPP $CPPFLAGS' 3571 ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5' 3572 ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' 3573 ac_compiler_gnu=$ac_cv_cxx_compiler_gnu 3574 3512 3575 3513 3576 # It seems that we need to cleanup something here for the Windows … … 5716 5779 *-*-irix6*) 5717 5780 # Find out which ABI we are using. 5718 echo '#line 57 18"configure"' > conftest.$ac_ext5781 echo '#line 5781 "configure"' > conftest.$ac_ext 5719 5782 if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5 5720 5783 (eval $ac_compile) 2>&5 … … 6850 6913 6851 6914 # Provide some information about the compiler. 6852 echo "$as_me:6 852:" \6915 echo "$as_me:6915:" \ 6853 6916 "checking for Fortran 77 compiler version" >&5 6854 6917 ac_compiler=`set X $ac_compile; echo $2` … … 7917 7980 -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ 7918 7981 -e 's:$: $lt_compiler_flag:'` 7919 (eval echo "\"\$as_me:79 19: $lt_compile\"" >&5)7982 (eval echo "\"\$as_me:7982: $lt_compile\"" >&5) 7920 7983 (eval "$lt_compile" 2>conftest.err) 7921 7984 ac_status=$? 7922 7985 cat conftest.err >&5 7923 echo "$as_me:79 23: \$? = $ac_status" >&57986 echo "$as_me:7986: \$? = $ac_status" >&5 7924 7987 if (exit $ac_status) && test -s "$ac_outfile"; then 7925 7988 # The compiler can only warn and ignore the option if not recognized … … 8185 8248 -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ 8186 8249 -e 's:$: $lt_compiler_flag:'` 8187 (eval echo "\"\$as_me:8 187: $lt_compile\"" >&5)8250 (eval echo "\"\$as_me:8250: $lt_compile\"" >&5) 8188 8251 (eval "$lt_compile" 2>conftest.err) 8189 8252 ac_status=$? 8190 8253 cat conftest.err >&5 8191 echo "$as_me:8 191: \$? = $ac_status" >&58254 echo "$as_me:8254: \$? = $ac_status" >&5 8192 8255 if (exit $ac_status) && test -s "$ac_outfile"; then 8193 8256 # The compiler can only warn and ignore the option if not recognized … … 8289 8352 -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ 8290 8353 -e 's:$: $lt_compiler_flag:'` 8291 (eval echo "\"\$as_me:8 291: $lt_compile\"" >&5)8354 (eval echo "\"\$as_me:8354: $lt_compile\"" >&5) 8292 8355 (eval "$lt_compile" 2>out/conftest.err) 8293 8356 ac_status=$? 8294 8357 cat out/conftest.err >&5 8295 echo "$as_me:8 295: \$? = $ac_status" >&58358 echo "$as_me:8358: \$? = $ac_status" >&5 8296 8359 if (exit $ac_status) && test -s out/conftest2.$ac_objext 8297 8360 then … … 10634 10697 lt_status=$lt_dlunknown 10635 10698 cat > conftest.$ac_ext <<EOF 10636 #line 106 36"configure"10699 #line 10699 "configure" 10637 10700 #include "confdefs.h" 10638 10701 … … 10734 10797 lt_status=$lt_dlunknown 10735 10798 cat > conftest.$ac_ext <<EOF 10736 #line 107 36"configure"10799 #line 10799 "configure" 10737 10800 #include "confdefs.h" 10738 10801 … … 13078 13141 -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ 13079 13142 -e 's:$: $lt_compiler_flag:'` 13080 (eval echo "\"\$as_me:13 080: $lt_compile\"" >&5)13143 (eval echo "\"\$as_me:13143: $lt_compile\"" >&5) 13081 13144 (eval "$lt_compile" 2>conftest.err) 13082 13145 ac_status=$? 13083 13146 cat conftest.err >&5 13084 echo "$as_me:13 084: \$? = $ac_status" >&513147 echo "$as_me:13147: \$? = $ac_status" >&5 13085 13148 if (exit $ac_status) && test -s "$ac_outfile"; then 13086 13149 # The compiler can only warn and ignore the option if not recognized … … 13182 13245 -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ 13183 13246 -e 's:$: $lt_compiler_flag:'` 13184 (eval echo "\"\$as_me:13 184: $lt_compile\"" >&5)13247 (eval echo "\"\$as_me:13247: $lt_compile\"" >&5) 13185 13248 (eval "$lt_compile" 2>out/conftest.err) 13186 13249 ac_status=$? 13187 13250 cat out/conftest.err >&5 13188 echo "$as_me:13 188: \$? = $ac_status" >&513251 echo "$as_me:13251: \$? = $ac_status" >&5 13189 13252 if (exit $ac_status) && test -s out/conftest2.$ac_objext 13190 13253 then … … 14752 14815 -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ 14753 14816 -e 's:$: $lt_compiler_flag:'` 14754 (eval echo "\"\$as_me:14 754: $lt_compile\"" >&5)14817 (eval echo "\"\$as_me:14817: $lt_compile\"" >&5) 14755 14818 (eval "$lt_compile" 2>conftest.err) 14756 14819 ac_status=$? 14757 14820 cat conftest.err >&5 14758 echo "$as_me:14 758: \$? = $ac_status" >&514821 echo "$as_me:14821: \$? = $ac_status" >&5 14759 14822 if (exit $ac_status) && test -s "$ac_outfile"; then 14760 14823 # The compiler can only warn and ignore the option if not recognized … … 14856 14919 -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ 14857 14920 -e 's:$: $lt_compiler_flag:'` 14858 (eval echo "\"\$as_me:14 858: $lt_compile\"" >&5)14921 (eval echo "\"\$as_me:14921: $lt_compile\"" >&5) 14859 14922 (eval "$lt_compile" 2>out/conftest.err) 14860 14923 ac_status=$? 14861 14924 cat out/conftest.err >&5 14862 echo "$as_me:14 862: \$? = $ac_status" >&514925 echo "$as_me:14925: \$? = $ac_status" >&5 14863 14926 if (exit $ac_status) && test -s out/conftest2.$ac_objext 14864 14927 then … … 17063 17126 -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ 17064 17127 -e 's:$: $lt_compiler_flag:'` 17065 (eval echo "\"\$as_me:17 065: $lt_compile\"" >&5)17128 (eval echo "\"\$as_me:17128: $lt_compile\"" >&5) 17066 17129 (eval "$lt_compile" 2>conftest.err) 17067 17130 ac_status=$? 17068 17131 cat conftest.err >&5 17069 echo "$as_me:17 069: \$? = $ac_status" >&517132 echo "$as_me:17132: \$? = $ac_status" >&5 17070 17133 if (exit $ac_status) && test -s "$ac_outfile"; then 17071 17134 # The compiler can only warn and ignore the option if not recognized … … 17331 17394 -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ 17332 17395 -e 's:$: $lt_compiler_flag:'` 17333 (eval echo "\"\$as_me:173 33: $lt_compile\"" >&5)17396 (eval echo "\"\$as_me:17396: $lt_compile\"" >&5) 17334 17397 (eval "$lt_compile" 2>conftest.err) 17335 17398 ac_status=$? 17336 17399 cat conftest.err >&5 17337 echo "$as_me:17 337: \$? = $ac_status" >&517400 echo "$as_me:17400: \$? = $ac_status" >&5 17338 17401 if (exit $ac_status) && test -s "$ac_outfile"; then 17339 17402 # The compiler can only warn and ignore the option if not recognized … … 17435 17498 -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ 17436 17499 -e 's:$: $lt_compiler_flag:'` 17437 (eval echo "\"\$as_me:17 437: $lt_compile\"" >&5)17500 (eval echo "\"\$as_me:17500: $lt_compile\"" >&5) 17438 17501 (eval "$lt_compile" 2>out/conftest.err) 17439 17502 ac_status=$? 17440 17503 cat out/conftest.err >&5 17441 echo "$as_me:17 441: \$? = $ac_status" >&517504 echo "$as_me:17504: \$? = $ac_status" >&5 17442 17505 if (exit $ac_status) && test -s out/conftest2.$ac_objext 17443 17506 then -
trunk/Clp/src/ClpSimplex.hpp
r1321 r1326 715 715 */ 716 716 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); 717 721 public: 718 722 /** For advanced use. When doing iterative solves things can get -
trunk/Clp/src/ClpSolve.cpp
r1325 r1326 3125 3125 ClpSimplex::solve(CoinStructuredModel * model) 3126 3126 { 3127 double time1 = CoinCpuTime();3128 3127 // analyze structure 3129 int numberColumns = model->numberColumns();3130 3128 int numberRowBlocks=model->numberRowBlocks(); 3131 3129 int numberColumnBlocks = model->numberColumnBlocks(); … … 3164 3162 } 3165 3163 } 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 } 3168 3183 bool masterColumns = (numberColumnBlocks==numberRowBlocks); 3169 3184 if ((masterColumns&&numberElementBlocks==2*numberRowBlocks-1) 3170 3185 ||(!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 3802 3222 } 3803 3223 /* This loads a model from a CoinStructuredModel object - returns number of errors. … … 4103 3523 return largest; 4104 3524 } 3525 // Solve using Dantzig-Wolfe decomposition and maybe in parallel 3526 int 3527 ClpSimplex::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 4177 int 4178 ClpSimplex::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.
