Changeset 383
 Timestamp:
 Jun 1, 2004 2:38:34 PM (15 years ago)
 Location:
 trunk
 Files:

 4 edited
Legend:
 Unmodified
 Added
 Removed

trunk/ClpSimplex.cpp
r372 r383 3487 3487 { 3488 3488 return ((ClpSimplexPrimalQuadratic *) this)>primalQuadratic(phase); 3489 } 3490 #include "ClpPredictorCorrector.hpp" 3491 #ifdef REAL_BARRIER 3492 #include "ClpCholeskyWssmp.hpp" 3493 #include "ClpCholeskyWssmpKKT.hpp" 3494 //#include "ClpCholeskyTaucs.hpp" 3495 #endif 3496 #include "ClpPresolve.hpp" 3497 /* Solves using barrier (assumes you have good cholesky factor code). 3498 Does crossover to simplex if asked*/ 3499 int 3500 ClpSimplex::barrier(bool crossover) 3501 { 3502 ClpSimplex * model2 = this; 3503 int savePerturbation=perturbation_; 3504 ClpInterior barrier; 3505 barrier.borrowModel(*model2); 3506 #ifdef REAL_BARRIER 3507 // uncomment this if you have Anshul Gupta's wsmp package 3508 ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(max(100,model2>numberRows()/10)); 3509 //ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(); 3510 //ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(max(100,model2>numberRows()/10)); 3511 // uncomment this if you have Sivan Toledo's Taucs package 3512 //ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs(); 3513 barrier.setCholesky(cholesky); 3514 #endif 3515 int numberRows = model2>numberRows(); 3516 int numberColumns = model2>numberColumns(); 3517 int saveMaxIts = model2>maximumIterations(); 3518 if (saveMaxIts<1000) { 3519 barrier.setMaximumBarrierIterations(saveMaxIts); 3520 model2>setMaximumIterations(1000000); 3521 } 3522 barrier.primalDual(); 3523 int barrierStatus=barrier.status(); 3524 double gap = barrier.complementarityGap(); 3525 // get which variables are fixed 3526 double * saveLower=NULL; 3527 double * saveUpper=NULL; 3528 ClpPresolve pinfo2; 3529 ClpSimplex * saveModel2=NULL; 3530 int numberFixed = barrier.numberFixed(); 3531 if (numberFixed*20>barrier.numberRows()&&numberFixed>5000&&crossover) { 3532 // may as well do presolve 3533 int numberRows = barrier.numberRows(); 3534 int numberColumns = barrier.numberColumns(); 3535 int numberTotal = numberRows+numberColumns; 3536 saveLower = new double [numberTotal]; 3537 saveUpper = new double [numberTotal]; 3538 memcpy(saveLower,barrier.columnLower(),numberColumns*sizeof(double)); 3539 memcpy(saveLower+numberColumns,barrier.rowLower(),numberRows*sizeof(double)); 3540 memcpy(saveUpper,barrier.columnUpper(),numberColumns*sizeof(double)); 3541 memcpy(saveUpper+numberColumns,barrier.rowUpper(),numberRows*sizeof(double)); 3542 barrier.fixFixed(); 3543 saveModel2=model2; 3544 } 3545 barrier.returnModel(*model2); 3546 double * rowPrimal = new double [numberRows]; 3547 double * columnPrimal = new double [numberColumns]; 3548 double * rowDual = new double [numberRows]; 3549 double * columnDual = new double [numberColumns]; 3550 // move solutions other way 3551 CoinMemcpyN(model2>primalRowSolution(), 3552 numberRows,rowPrimal); 3553 CoinMemcpyN(model2>dualRowSolution(), 3554 numberRows,rowDual); 3555 CoinMemcpyN(model2>primalColumnSolution(), 3556 numberColumns,columnPrimal); 3557 CoinMemcpyN(model2>dualColumnSolution(), 3558 numberColumns,columnDual); 3559 if (saveModel2) { 3560 // do presolve 3561 model2 = pinfo2.presolvedModel(*model2,1.0e8, 3562 false,5,true); 3563 } 3564 if (barrierStatus<4&&crossover) { 3565 // make sure no status left 3566 model2>createStatus(); 3567 // solve 3568 model2>setPerturbation(100); 3569 // throw some into basis 3570 { 3571 int numberRows = model2>numberRows(); 3572 int numberColumns = model2>numberColumns(); 3573 double * dsort = new double[numberColumns]; 3574 int * sort = new int[numberColumns]; 3575 int n=0; 3576 const double * columnLower = model2>columnLower(); 3577 const double * columnUpper = model2>columnUpper(); 3578 const double * primalSolution = model2>primalColumnSolution(); 3579 double tolerance = 10.0*primalTolerance_; 3580 int i; 3581 for ( i=0;i<numberRows;i++) 3582 model2>setRowStatus(i,superBasic); 3583 for ( i=0;i<numberColumns;i++) { 3584 double distance = min(columnUpper[i]primalSolution[i], 3585 primalSolution[i]columnLower[i]); 3586 if (distance>tolerance) { 3587 dsort[n]=distance; 3588 sort[n++]=i; 3589 model2>setStatus(i,superBasic); 3590 } else if (distance>primalTolerance_) { 3591 model2>setStatus(i,superBasic); 3592 } else if (primalSolution[i]<=columnLower[i]+primalTolerance_) { 3593 model2>setStatus(i,atLowerBound); 3594 } else { 3595 model2>setStatus(i,atUpperBound); 3596 } 3597 } 3598 CoinSort_2(dsort,dsort+n,sort); 3599 n = min(numberRows,n); 3600 for ( i=0;i<n;i++) { 3601 int iColumn = sort[i]; 3602 model2>setStatus(iColumn,basic); 3603 } 3604 delete [] sort; 3605 delete [] dsort; 3606 } 3607 if (gap<1.0e3*((double) (numberRows+numberColumns))) { 3608 int numberRows = model2>numberRows(); 3609 int numberColumns = model2>numberColumns(); 3610 // just primal values pass 3611 model2>primal(2); 3612 // save primal solution and copy back dual 3613 CoinMemcpyN(model2>primalRowSolution(), 3614 numberRows,rowPrimal); 3615 CoinMemcpyN(rowDual, 3616 numberRows,model2>dualRowSolution()); 3617 CoinMemcpyN(model2>primalColumnSolution(), 3618 numberColumns,columnPrimal); 3619 CoinMemcpyN(columnDual, 3620 numberColumns,model2>dualColumnSolution()); 3621 //model2>primal(1); 3622 // clean up reduced costs and flag variables 3623 { 3624 double * dj = model2>dualColumnSolution(); 3625 double * cost = model2>objective(); 3626 double * saveCost = new double[numberColumns]; 3627 memcpy(saveCost,cost,numberColumns*sizeof(double)); 3628 double * saveLower = new double[numberColumns]; 3629 double * lower = model2>columnLower(); 3630 memcpy(saveLower,lower,numberColumns*sizeof(double)); 3631 double * saveUpper = new double[numberColumns]; 3632 double * upper = model2>columnUpper(); 3633 memcpy(saveUpper,upper,numberColumns*sizeof(double)); 3634 int i; 3635 double tolerance = 10.0*dualTolerance_; 3636 for ( i=0;i<numberColumns;i++) { 3637 if (model2>getStatus(i)==basic) { 3638 dj[i]=0.0; 3639 } else if (model2>getStatus(i)==atLowerBound) { 3640 if (optimizationDirection_*dj[i]<tolerance) { 3641 if (optimizationDirection_*dj[i]<0.0) { 3642 //if (dj[i]<1.0e3) 3643 //printf("bad dj at lb %d %g\n",i,dj[i]); 3644 cost[i] = dj[i]; 3645 dj[i]=0.0; 3646 } 3647 } else { 3648 upper[i]=lower[i]; 3649 } 3650 } else if (model2>getStatus(i)==atUpperBound) { 3651 if (optimizationDirection_*dj[i]>tolerance) { 3652 if (optimizationDirection_*dj[i]>0.0) { 3653 //if (dj[i]>1.0e3) 3654 //printf("bad dj at ub %d %g\n",i,dj[i]); 3655 cost[i] = dj[i]; 3656 dj[i]=0.0; 3657 } 3658 } else { 3659 lower[i]=upper[i]; 3660 } 3661 } 3662 } 3663 // just dual values pass 3664 //model2>setLogLevel(63); 3665 //model2>setFactorizationFrequency(1); 3666 model2>dual(2); 3667 memcpy(cost,saveCost,numberColumns*sizeof(double)); 3668 delete [] saveCost; 3669 memcpy(lower,saveLower,numberColumns*sizeof(double)); 3670 delete [] saveLower; 3671 memcpy(upper,saveUpper,numberColumns*sizeof(double)); 3672 delete [] saveUpper; 3673 } 3674 // and finish 3675 // move solutions 3676 CoinMemcpyN(rowPrimal, 3677 numberRows,model2>primalRowSolution()); 3678 CoinMemcpyN(columnPrimal, 3679 numberColumns,model2>primalColumnSolution()); 3680 } 3681 model2>primal(1); 3682 } else if (barrierStatus==4&&crossover) { 3683 // memory problems 3684 model2>setPerturbation(savePerturbation); 3685 model2>createStatus(); 3686 model2>dual(); 3687 } 3688 model2>setMaximumIterations(saveMaxIts); 3689 delete [] rowPrimal; 3690 delete [] columnPrimal; 3691 delete [] rowDual; 3692 delete [] columnDual; 3693 if (saveLower) { 3694 pinfo2.postsolve(true); 3695 delete model2; 3696 model2=saveModel2; 3697 int numberRows = model2>numberRows(); 3698 int numberColumns = model2>numberColumns(); 3699 memcpy(model2>columnLower(),saveLower,numberColumns*sizeof(double)); 3700 memcpy(model2>rowLower(),saveLower+numberColumns,numberRows*sizeof(double)); 3701 delete [] saveLower; 3702 memcpy(model2>columnUpper(),saveUpper,numberColumns*sizeof(double)); 3703 memcpy(model2>rowUpper(),saveUpper+numberColumns,numberRows*sizeof(double)); 3704 delete [] saveUpper; 3705 model2>primal(1); 3706 } 3707 model2>setPerturbation(savePerturbation); 3708 return model2>status(); 3489 3709 } 3490 3710 /* For strong branching. On input lower and upper are new bounds 
trunk/Test/ClpMain.cpp
r374 r383 14 14 #include "CoinPragma.hpp" 15 15 #include "CoinHelperFunctions.hpp" 16 #define CLPVERSION "0.99. 6"16 #define CLPVERSION "0.99.7" 17 17 18 18 #include "CoinMpsIO.hpp" … … 50 50 #endif 51 51 52 int mainTest (int argc, const char *argv[], bool doDual,52 int mainTest (int argc, const char *argv[],int algorithm, 53 53 ClpSimplex empty, bool doPresolve,int doIdiot); 54 54 enum ClpParameterType { 55 GENERALQUERY=100, 55 GENERALQUERY=100,FULLGENERALQUERY, 56 56 57 57 DUALTOLERANCE=1,PRIMALTOLERANCE,DUALBOUND,PRIMALWEIGHT,MAXTIME,OBJSCALE, … … 66 66 DIRECTORY=301,IMPORT,EXPORT,RESTORE,SAVE,DUALSIMPLEX,PRIMALSIMPLEX, 67 67 MAXIMIZE,MINIMIZE,EXIT,STDIN,UNITTEST,NETLIB_DUAL,NETLIB_PRIMAL,SOLUTION, 68 TIGHTEN,FAKEBOUND,HELP,PLUSMINUS,NETWORK,ALLSLACK,REVERSE,BARRIER, 68 TIGHTEN,FAKEBOUND,HELP,PLUSMINUS,NETWORK,ALLSLACK,REVERSE,BARRIER,NETLIB_BARRIER, 69 69 70 70 INVALID=1000 … … 897 897 ClpItem("?","For help",GENERALQUERY,1,false); 898 898 parameters[numberParameters++]= 899 ClpItem("???","For help",FULLGENERALQUERY,1,false); 900 parameters[numberParameters++]= 899 901 ClpItem("","From stdin", 900 902 STDIN,299,false); … … 922 924 ( 923 925 "This command solves the current model using the primal dual predictor \ 924 corrector algorithm. This is not a sophisticated version just something JJF\ 925 knocked up\ 926 ** another slight drawback (early December) is that it does not work" 926 corrector algorithm. This is not a sophisticated version just something JJF \ 927 knocked up" 927 928 928 929 ); … … 1150 1151 The user can set options before e.g. clp presolve off netlib" 1151 1152 ); 1153 #ifdef REAL_BARRIER 1154 parameters[numberParameters++]= 1155 ClpItem("netlibB!arrier","Solve entire netlib test set with barrier", 1156 NETLIB_BARRIER); 1157 parameters[numberParameters1].setLonghelp 1158 ( 1159 "This exercises the unit test for clp and then solves the netlib test set using barrier.\ 1160 The user can set options before e.g. clp presolve off netlib" 1161 ); 1162 #endif 1152 1163 parameters[numberParameters++]= 1153 1164 ClpItem("netlibP!rimal","Solve entire netlib test set (primal)", … … 1444 1455 // see if ? at end 1445 1456 int numberQuery=0; 1446 if (field!="?" ) {1457 if (field!="?"&&field!="???") { 1447 1458 int length = field.length(); 1448 1459 int i; … … 1500 1511 int type = parameters[iParam].indexNumber(); 1501 1512 if (parameters[iParam].displayThis()&&type>=limits[iType] 1513 &&type<limits[iType+1]) { 1514 if (!across) 1515 std::cout<<" "; 1516 std::cout<<parameters[iParam].matchName()<<" "; 1517 across++; 1518 if (across==maxAcross) { 1519 std::cout<<std::endl; 1520 across=0; 1521 } 1522 } 1523 } 1524 if (across) 1525 std::cout<<std::endl; 1526 } 1527 } else if (type==FULLGENERALQUERY) { 1528 std::cout<<"Full list of ommands is:"<<std::endl; 1529 int maxAcross=5; 1530 int limits[]={1,101,201,301,401}; 1531 std::vector<std::string> types; 1532 types.push_back("Double parameters:"); 1533 types.push_back("Int parameters:"); 1534 types.push_back("Keyword parameters and others:"); 1535 types.push_back("Actions:"); 1536 int iType; 1537 for (iType=0;iType<4;iType++) { 1538 int across=0; 1539 std::cout<<types[iType]<<std::endl; 1540 for ( iParam=0; iParam<numberParameters; iParam++ ) { 1541 int type = parameters[iParam].indexNumber(); 1542 if (type>=limits[iType] 1502 1543 &&type<limits[iType+1]) { 1503 1544 if (!across) … … 2130 2171 break; 2131 2172 case NETLIB_DUAL: 2173 case NETLIB_BARRIER: 2132 2174 case NETLIB_PRIMAL: 2133 2175 { … … 2142 2184 nFields=4; 2143 2185 } 2144 if (type==NETLIB_DUAL) 2186 int algorithm; 2187 if (type==NETLIB_DUAL) { 2145 2188 std::cerr<<"Doing netlib with dual agorithm"<<std::endl; 2146 else 2189 algorithm =0; 2190 } else if (type==NETLIB_BARRIER) { 2191 std::cerr<<"Doing netlib with barrier agorithm"<<std::endl; 2192 algorithm =2; 2193 } else { 2147 2194 std::cerr<<"Doing netlib with primal agorithm"<<std::endl; 2148 mainTest(nFields,fields,(type==NETLIB_DUAL),models[iModel], 2195 algorithm=1; 2196 } 2197 mainTest(nFields,fields,algorithm,models[iModel], 2149 2198 (preSolve!=0),doIdiot); 2150 2199 } 
trunk/Test/unitTest.cpp
r345 r383 56 56 // All parameters are optional. 57 57 // 58 int mainTest (int argc, const char *argv[], bool doDual,58 int mainTest (int argc, const char *argv[],int algorithm, 59 59 ClpSimplex empty, bool doPresolve,int doIdiot) 60 60 { … … 275 275 model2>numberRows()>5000) 276 276 model2>factorization()>maximumPivots(100+model2>numberRows()/100); 277 if ( doDual) {277 if (algorithm==0) { 278 278 // faster if bounds tightened 279 279 int numberInfeasibilities = model2>tightenPrimalBounds(); … … 284 284 //model2>crash(1000,1); 285 285 model2>dual(); 286 } else if (algorithm==2) { 287 model2>barrier(); 286 288 } else { 287 289 if (doIdiot>0) { … … 322 324 } 323 325 } else { 324 if ( doDual) {326 if (algorithm==0) { 325 327 if (doIdiot<0) 326 328 solution.crash(1000,1); 327 329 solution.dual(); 330 } else if (algorithm==2) { 331 solution.barrier(); 328 332 } else { 329 333 if (doIdiot>0) { 
trunk/include/ClpSimplex.hpp
r372 r383 178 178 /// Solves quadratic using Dantzig's algorithm  primal 179 179 int quadraticPrimal(int phase=2); 180 /** Solves using barrier (assumes you have good cholesky factor code). 181 Does crossover to simplex if asked*/ 182 int barrier(bool crossover=true); 180 183 /** Dual ranging. 181 184 This computes increase/decrease in cost for each given variable and corresponding
Note: See TracChangeset
for help on using the changeset viewer.