Changeset 378


Ignore:
Timestamp:
May 27, 2004 5:10:30 PM (14 years ago)
Author:
forrest
Message:

Interior

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpInterior.cpp

    r375 r378  
    640640  assert (!dj_);
    641641  dj_ = new double [nTotal];
    642   delete [] status_;
    643   status_ = new unsigned char [numberRows_+numberColumns_];
     642  if (!status_)
     643    status_ = new unsigned char [numberRows_+numberColumns_];
    644644  memset(status_,0,numberRows_+numberColumns_);
    645645  return goodMatrix;
     
    943943  sumPrimalInfeasibilities_=0.0;
    944944  sumDualInfeasibilities_=0.0;
    945   double dualTolerance =  dblParam_[ClpDualTolerance];
     945  double dualTolerance =  10.0*dblParam_[ClpDualTolerance];
    946946  double primalTolerance =  dblParam_[ClpPrimalTolerance];
    947947  worstComplementarity_=0.0;
     
    10271027  cholesky_= cholesky;
    10281028}
     1029/* Borrow model.  This is so we dont have to copy large amounts
     1030   of data around.  It assumes a derived class wants to overwrite
     1031   an empty model with a real one - while it does an algorithm.
     1032   This is same as ClpModel one. */
     1033void
     1034ClpInterior::borrowModel(ClpModel & otherModel)
     1035{
     1036  ClpModel::borrowModel(otherModel);
     1037}
     1038/* Return model - updates any scalars */
     1039void
     1040ClpInterior::returnModel(ClpModel & otherModel)
     1041{
     1042  ClpModel::returnModel(otherModel);
     1043}
  • trunk/ClpPredictorCorrector.cpp

    r375 r378  
    3131static double eDiagonalCaution=1.0e18;
    3232static double eExtra=1.0e-12;
     33static double eFree =1.0e3;
    3334
    3435// main function
     
    412413    //save information
    413414    double product=affineProduct();
     415    //#define ALWAYS
     416#ifndef ALWAYS
    414417#if 0
    415418    //can we do corrector step?
     
    461464      phase=1;
    462465    }
     466#else
     467    phase=1;
     468#endif
    463469    if (goodMove&&doCorrector) {
    464470      //set up for next step
     
    489495        findStepLength(phase);
    490496        nextGap = complementarityGap(nextNumber,nextNumberItems,1);
     497#ifndef ALWAYS
    491498        if (numberIterations_>=-77) {
    492499          goodMove=checkGoodMove(true,bestNextGap);
     
    500507              // Back to affine
    501508              phase=0;
    502               // Try primal dual step instead - but with small mu
    503               phase=2;
    504               double floatNumber;
    505               floatNumber = 2.0*numberComplementarityPairs_;
    506               mu_=complementarityGap_/floatNumber;
    507               double mu1=mu_;
    508               double phi;
    509               if (numberComplementarityPairs_<=500) {
    510                 phi=pow((double) numberComplementarityPairs_,2.0);
    511               } else {
    512                 phi=pow((double) numberComplementarityPairs_,1.5);
    513                 if (phi<500.0*500.0) {
    514                   phi=500.0*500.0;
    515                 }
     509              if (numberIterations_<100) {
     510                // Try primal dual step instead - but with small mu
     511                phase=2;
     512                double floatNumber;
     513                floatNumber = 2.0*numberComplementarityPairs_;
     514                mu_=complementarityGap_/floatNumber;
     515                double mu1=mu_;
     516                double phi;
     517                if (numberComplementarityPairs_<=500) {
     518                  phi=pow((double) numberComplementarityPairs_,2.0);
     519                } else {
     520                  phi=pow((double) numberComplementarityPairs_,1.5);
     521                  if (phi<500.0*500.0) {
     522                    phi=500.0*500.0;
     523                  }
     524                }
     525                mu_=complementarityGap_/phi;
     526                //printf("pd mu %g, alternate %g, smallest %g\n",
     527                //     mu_,mu1,smallestPrimalDualMu);
     528                mu_ = sqrt(mu_*mu1);
     529                mu_=mu1*0.8;
    516530              }
    517               mu_=complementarityGap_/phi;
    518               //printf("pd mu %g, alternate %g, smallest %g\n",
    519               //     mu_,mu1,smallestPrimalDualMu);
    520               mu_ = sqrt(mu_*mu1);
    521               mu_=mu1*0.8;
    522531              setupForSolve(phase);
    523532              directionAccuracy=findDirectionVector(phase);
     
    530539          goodMove=true;
    531540        }
     541#endif
    532542      }
    533543    }
     
    11581168    if (!flagged(iColumn)) {
    11591169      double deltaX = deltaX_[iColumn];
    1160       if (lowerBound(iColumn)) {
    1161         double deltaSL = rhsL_[iColumn]+deltaX;
    1162         double slack = lowerSlack[iColumn]+extra;
    1163         deltaSL_[iColumn]=deltaSL;
    1164         deltaZ_[iColumn]=(rhsZ_[iColumn]-zVec[iColumn]*deltaSL)/slack;
    1165       }
    1166       if (upperBound(iColumn)) {
    1167         double deltaSU = rhsU_[iColumn]-deltaX;
    1168         double slack = upperSlack[iColumn]+extra;
    1169         deltaSU_[iColumn]=deltaSU;
    1170         deltaT_[iColumn]=(rhsT_[iColumn]-tVec[iColumn]*deltaSU)/slack;
     1170      if (lowerBound(iColumn)||upperBound(iColumn)) {
     1171        if (lowerBound(iColumn)) {
     1172          double deltaSL = rhsL_[iColumn]+deltaX;
     1173          double slack = lowerSlack[iColumn]+extra;
     1174          deltaSL_[iColumn]=deltaSL;
     1175          deltaZ_[iColumn]=(rhsZ_[iColumn]-zVec[iColumn]*deltaSL)/slack;
     1176        }
     1177        if (upperBound(iColumn)) {
     1178          double deltaSU = rhsU_[iColumn]-deltaX;
     1179          double slack = upperSlack[iColumn]+extra;
     1180          deltaSU_[iColumn]=deltaSU;
     1181          deltaT_[iColumn]=(rhsT_[iColumn]-tVec[iColumn]*deltaSU)/slack;
     1182        }
     1183      } else {
     1184        // free
     1185        double gap=fabs(solution_[iColumn]);
     1186        double multiplier=1.0/gap;
     1187        if (gap<1.0) {
     1188          multiplier=1.0;
     1189        }
     1190        deltaZ_[iColumn] = -zVec[iColumn] - multiplier*deltaX*zVec[iColumn];
     1191        deltaT_[iColumn] = -tVec[iColumn] + multiplier*deltaX*tVec[iColumn];
    11711192      }
    11721193    }
     
    15311552  for ( iColumn=0;iColumn<numberTotal;iColumn++) {
    15321553    if (!flagged(iColumn)) {
    1533       double newValue=solution_[iColumn];
    1534       double lowerValue=lower_[iColumn];
    1535       double upperValue=upper_[iColumn];
    1536       if (newValue>lowerValue+largeGap&&newValue<upperValue-largeGap) {
     1554      if (!lowerBound(iColumn)&&!upperBound(iColumn)) {
    15371555        clearFixedOrFree(iColumn);
    15381556        setLowerBound(iColumn);
    15391557        setUpperBound(iColumn);
    15401558        double objectiveValue = cost_[iColumn];
    1541         lowerValue=max(lowerValue,newValue-largeGap);
    1542         upperValue=min(upperValue,newValue+largeGap);
     1559        double newValue=solution_[iColumn];
     1560        double lowerValue=max(lower_[iColumn],newValue-largeGap);
     1561        double upperValue=min(upper_[iColumn],newValue+largeGap);
    15431562        lower_[iColumn]=lowerValue;
    15441563        upper_[iColumn]=upperValue;
     
    16551674        } else {
    16561675          double change;
    1657           change =primal[iColumn]+deltaX[iColumn]-lowerSlack[iColumn]-lower[iColumn];
     1676          if (!fakeLower(iColumn)) {
     1677            change =primal[iColumn]+deltaX[iColumn]-lowerSlack[iColumn]-lower[iColumn];
     1678          } else {
     1679            change =deltaX[iColumn];
     1680          }
    16581681          dualValue=zVec[iColumn]+actualDualStep_*deltaZ[iColumn];
    16591682          primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
     
    16891712        } else {
    16901713          double change;
    1691           change =upper[iColumn]-primal[iColumn]-deltaX[iColumn]-upperSlack[iColumn];
     1714          if (!fakeUpper(iColumn)) {
     1715            change =upper[iColumn]-primal[iColumn]-deltaX[iColumn]-upperSlack[iColumn];
     1716          } else {
     1717            change =-deltaX[iColumn];
     1718          }
    16921719          dualValue=tVec[iColumn]+actualDualStep_*deltaT[iColumn];
    16931720          primalValue=upperSlack[iColumn]+actualPrimalStep_*change;
     
    17601787        rhsC_[iColumn] = dj[iColumn]-zVec[iColumn]+tVec[iColumn];
    17611788        if (lowerBound(iColumn)) {
    1762           rhsZ_[iColumn] = -zVec[iColumn]*(lowerSlack[iColumn]+extra);
    1763           rhsL_[iColumn] = min(0.0,primal[iColumn]-lowerSlack[iColumn]-lower[iColumn]);
     1789          if (!fakeLower(iColumn)) {
     1790            rhsZ_[iColumn] = -zVec[iColumn]*(lowerSlack[iColumn]+extra);
     1791            rhsL_[iColumn] = min(0.0,primal[iColumn]-lowerSlack[iColumn]-lower[iColumn]);
     1792          }
     1793        } else if (upperBound(iColumn)){
     1794          assert (!zVec[iColumn]);
    17641795        }
    17651796        if (upperBound(iColumn)) {
    1766           rhsT_[iColumn] = -tVec[iColumn]*(upperSlack[iColumn]+extra);
    1767           rhsU_[iColumn] = min(0.0,-(primal[iColumn] + upperSlack[iColumn] - upper[iColumn]));
     1797          if (!fakeUpper(iColumn)) {
     1798            rhsT_[iColumn] = -tVec[iColumn]*(upperSlack[iColumn]+extra);
     1799            rhsU_[iColumn] = min(0.0,-(primal[iColumn] + upperSlack[iColumn] - upper[iColumn]));
     1800          }
     1801        } else if (lowerBound(iColumn)) {
     1802          assert (!tVec[iColumn]);
    17681803        }
    17691804      }
     
    18281863          // To bring in line with OSL
    18291864          rhsZ_[iColumn] -= deltaZ_[iColumn]*rhsL_[iColumn];
     1865          //rhsZ_[iColumn] = mu_ -deltaZ_[iColumn]*rhsL_[iColumn]
     1866          //- deltaZ_[iColumn]*deltaX_[iColumn];
    18301867        }
    18311868        if (upperBound(iColumn)) {
     
    18341871          // To bring in line with OSL
    18351872          rhsT_[iColumn] -= deltaT_[iColumn]*rhsU_[iColumn];
     1873          //rhsT_[iColumn] = mu_ -deltaT_[iColumn]*rhsU_[iColumn]
     1874          //+deltaT_[iColumn]*deltaX_[iColumn];
    18361875        }
    18371876      }
     
    18761915    for (iColumn=0;iColumn<numberTotal;iColumn++) {
    18771916      rhsC_[iColumn]=0.0;
     1917      //rhsU_[iColumn]=0.0;
     1918      //rhsL_[iColumn]=0.0;
    18781919      rhsZ_[iColumn]=0.0;
    18791920      rhsT_[iColumn]=0.0;
     
    18811922        rhsC_[iColumn] = dj[iColumn]-zVec[iColumn]+tVec[iColumn];
    18821923        if (lowerBound(iColumn)) {
    1883           rhsZ_[iColumn] = mu_ - zVec[iColumn]*(lowerSlack[iColumn]+extra);
     1924          if (!fakeLower(iColumn)) {
     1925            rhsZ_[iColumn] = mu_ - zVec[iColumn]*(lowerSlack[iColumn]+extra);
     1926            //rhsL_[iColumn] = primal[iColumn]-lowerSlack[iColumn]-lower[iColumn];
     1927          }
    18841928        }
    18851929        if (upperBound(iColumn)) {
    1886           rhsT_[iColumn] = mu_ - tVec[iColumn]*(upperSlack[iColumn]+extra);
     1930          if (!fakeUpper(iColumn)) {
     1931            rhsT_[iColumn] = mu_ - tVec[iColumn]*(upperSlack[iColumn]+extra);
     1932            //rhsU_[iColumn] = -(primal[iColumn] + upperSlack[iColumn] - upper[iColumn]);
     1933          }
    18871934        }
    18881935      }
     
    19251972              assert (value<0.0);
    19261973            }
     1974            //#define AGAIN
     1975#ifdef AGAIN
     1976            //rhsZ_[iColumn] = mu_ -zVec[iColumn]*(lowerSlack[iColumn]+extra)
     1977            //- deltaZ_[iColumn]*deltaX_[iColumn];
     1978            // To bring in line with OSL
     1979            rhsZ_[iColumn] += deltaZ_[iColumn]*rhsL_[iColumn];
     1980#endif
    19271981            rhsZ_[iColumn] += value;
    19281982          } 
     
    19482002              assert (value<0.0);
    19492003            }
     2004#ifdef AGAIN
     2005            //rhsT_[iColumn] = mu_ -tVec[iColumn]*(upperSlack[iColumn]+extra)
     2006            //+deltaT_[iColumn]*deltaX_[iColumn];
     2007            // To bring in line with OSL
     2008            rhsT_[iColumn] -= deltaT_[iColumn]*rhsU_[iColumn];
     2009#endif
    19502010            rhsT_[iColumn] += value;
    19512011          }
     
    22882348  }
    22892349  //qDiagonal *= 1.0e2;
     2350  double widenGap=1.0e1;
    22902351  //largest allowable ratio of lowerSlack/zVec (etc)
    22912352  double largestRatio;
     
    23542415        largerzw=tValue;
    23552416      }
    2356       bool fakeOldBounds=false;
    2357       bool fakeNewBounds=false;
    23582417      double trueLower;
    23592418      double trueUpper;
     
    23652424        trueUpper = rowUpper_[iColumn-numberColumns_];
    23662425      }
    2367       if (oldPrimal>trueLower+largeGap2&&
    2368           oldPrimal<trueUpper-largeGap2)
    2369         fakeOldBounds=true;
    2370       if (newPrimal>trueLower+largeGap2&&
    2371           newPrimal<trueUpper-largeGap2)
    2372         fakeNewBounds=true;
    2373       if (fakeOldBounds) {
    2374         if (fakeNewBounds) {
    2375           lower_[iColumn]=newPrimal-largeGap2;
    2376           lowerSlack[iColumn] = largeGap2;
    2377           upper_[iColumn]=newPrimal+largeGap2;
    2378           upperSlack[iColumn] = largeGap2;
    2379         } else {
    2380           lower_[iColumn]=trueLower;
    2381           setLowerBound(iColumn);
    2382           lowerSlack[iColumn] = max(newPrimal-trueLower,1.0);
    2383           upper_[iColumn]=trueUpper;
    2384           setUpperBound(iColumn);
    2385           upperSlack[iColumn] = max(trueUpper-newPrimal,1.0);
    2386         }
    2387       } else if (fakeNewBounds) {
    2388         lower_[iColumn]=newPrimal-largeGap2;
    2389         lowerSlack[iColumn] = largeGap2;
    2390         upper_[iColumn]=newPrimal+largeGap2;
    2391         upperSlack[iColumn] = largeGap2;
    2392         // so we can just have one test
    2393         fakeOldBounds=true;
    2394       }
    23952426      if (lowerBound(iColumn)) {
    23962427        double oldSlack = lowerSlack[iColumn];
    23972428        double newSlack;
    2398         newSlack=
    2399           lowerSlack[iColumn]+actualPrimalStep_*(oldPrimal-oldSlack
    2400                                                  + thisWeight-lower[iColumn]);
    2401         if (fakeOldBounds)
    2402           newSlack = lowerSlack[iColumn];
     2429        if (!fakeLower(iColumn)) {
     2430          newSlack=
     2431                   lowerSlack[iColumn]+actualPrimalStep_*(oldPrimal-oldSlack
     2432                + thisWeight-lower[iColumn]);
     2433        } else {
     2434          newSlack= lowerSlack[iColumn]+actualPrimalStep_*thisWeight;
     2435          if (newSlack<0.0) {
     2436            abort();
     2437          }
     2438        }
     2439        if (trueLower!=lower_[iColumn]) {
     2440          if (newPrimal>trueLower+largeGap2) {
     2441            lower_[iColumn]=newPrimal-largeGap2;
     2442            newSlack = largeGap2;
     2443          } else {
     2444            lower_[iColumn]=trueLower;
     2445            newSlack = max(0.5*newSlack,newPrimal-trueLower);
     2446          }
     2447        }
    24032448        double epsilon = fabs(newSlack)*epsilonBase;
    24042449        if (epsilon>1.0e-5) {
     
    24232468            smallGap2=0.1*fabs(newPrimal);
    24242469          }
    2425           double larger;
    2426           if (newSlack>feasibleSlack) {
    2427             larger=newSlack;
    2428           } else {
    2429             larger=feasibleSlack;
    2430           }
    2431           if (fabs(feasibleSlack-newSlack)<1.0e-6*larger) {
    2432             newSlack=feasibleSlack;
    2433           }
     2470          if (!fakeLower(iColumn)) {
     2471            double larger;
     2472            if (newSlack>feasibleSlack) {
     2473              larger=newSlack;
     2474            } else {
     2475              larger=feasibleSlack;
     2476            }
     2477            if (fabs(feasibleSlack-newSlack)<1.0e-6*larger) {
     2478              newSlack=feasibleSlack;
     2479            }
     2480            //set FAKE here
     2481            if (newSlack>zValue2*largestRatio&&newSlack>smallGap2) {
     2482              setFakeLower(iColumn);
     2483              newSlack=zValue2*largestRatio;
     2484              if (newSlack<smallGap2) {
     2485                newSlack=smallGap2;
     2486              }
     2487              numberDecreased++;
     2488            }
     2489          } else {
     2490            newSlack=zValue2*largestRatio;
     2491            if (newSlack<smallGap2) {
     2492              newSlack=smallGap2;
     2493            }
     2494            if (newSlack>largeGap) {
     2495              //increase up to smaller of z.. and largeGap
     2496              newSlack=largeGap;
     2497            }
     2498            if (newSlack>widenGap*oldSlack) {
     2499              newSlack=widenGap*oldSlack;
     2500              numberIncreased++;
     2501              //cout<<"wider "<<j<<" "<<newSlack<<" "<<oldSlack<<" "<<feasibleSlack<<endl;
     2502            }
     2503            if (newSlack>feasibleSlack) {
     2504              newSlack=feasibleSlack;
     2505              clearFakeLower(iColumn);
     2506            }
     2507          }
    24342508        }
    24352509        if (zVec[iColumn]>dualTolerance) {
     
    24402514          smallerSlack=newSlack;
    24412515        }
    2442         double infeasibility = fabs(newPrimal-lowerSlack[iColumn]-lower[iColumn]);
    2443         if (infeasibility>maximumBoundInfeasibility) {
    2444           maximumBoundInfeasibility=infeasibility;
    2445         }
     2516        if (!fakeLower(iColumn)) {
     2517          double infeasibility = fabs(newPrimal-lowerSlack[iColumn]-lower[iColumn]);
     2518          if (infeasibility>maximumBoundInfeasibility) {
     2519            maximumBoundInfeasibility=infeasibility;
     2520          }
     2521        }
    24462522        if (lowerSlack[iColumn]<=kill&&fabs(newPrimal-lower[iColumn])<=kill) {
    24472523          //may be better to leave at value?
     
    24592535        double oldSlack = upperSlack[iColumn];
    24602536        double newSlack;
    2461         newSlack=
    2462           upperSlack[iColumn]+actualPrimalStep_*(-oldPrimal-oldSlack
    2463                                                  - thisWeight+upper[iColumn]);
    2464         if (fakeOldBounds)
    2465           newSlack = upperSlack[iColumn];
     2537        if (!fakeUpper(iColumn)) {
     2538          newSlack=
     2539                   upperSlack[iColumn]+actualPrimalStep_*(-oldPrimal-oldSlack
     2540                - thisWeight+upper[iColumn]);
     2541        } else {
     2542          newSlack= upperSlack[iColumn]-actualPrimalStep_*thisWeight;
     2543        }
     2544        if (trueUpper!=upper_[iColumn]) {
     2545          if (newPrimal<trueUpper-largeGap2) {
     2546            upper_[iColumn]=newPrimal+largeGap2;
     2547            newSlack = largeGap2;
     2548          } else {
     2549            upper_[iColumn]=trueUpper;
     2550            newSlack = max(0.5*newSlack,trueUpper-newPrimal);
     2551          }
     2552        }
    24662553        double epsilon = fabs(newSlack)*epsilonBase;
    24672554        if (epsilon>1.0e-5) {
     
    24862573            smallGap2=0.1*fabs(newPrimal);
    24872574          }
    2488           double larger;
    2489           if (newSlack>feasibleSlack) {
    2490             larger=newSlack;
    2491           } else {
    2492             larger=feasibleSlack;
    2493           }
    2494           if (fabs(feasibleSlack-newSlack)<1.0e-6*larger) {
    2495             newSlack=feasibleSlack;
    2496           }
     2575          if (!fakeUpper(iColumn)) {
     2576            double larger;
     2577            if (newSlack>feasibleSlack) {
     2578              larger=newSlack;
     2579            } else {
     2580              larger=feasibleSlack;
     2581            }
     2582            if (fabs(feasibleSlack-newSlack)<1.0e-6*larger) {
     2583              newSlack=feasibleSlack;
     2584            }
     2585            //set FAKE here
     2586            if (newSlack>tValue2*largestRatio&&newSlack>smallGap2) {
     2587              setFakeUpper(iColumn);
     2588              newSlack=tValue2*largestRatio;
     2589              if (newSlack<smallGap2) {
     2590                newSlack=smallGap2;
     2591              }
     2592              numberDecreased++;
     2593            }
     2594          } else {
     2595            newSlack=tValue2*largestRatio;
     2596            if (newSlack<smallGap2) {
     2597              newSlack=smallGap2;
     2598            }
     2599            if (newSlack>largeGap) {
     2600              //increase up to smaller of w.. and largeGap
     2601              newSlack=largeGap;
     2602            }
     2603            if (newSlack>widenGap*oldSlack) {
     2604              numberIncreased++;
     2605              //cout<<"wider "<<-j<<" "<<newSlack<<" "<<oldSlack<<" "<<feasibleSlack<<endl;
     2606            }
     2607            if (newSlack>feasibleSlack) {
     2608              newSlack=feasibleSlack;
     2609              clearFakeUpper(iColumn);
     2610            }
     2611          }
    24972612        }
    24982613        if (tVec[iColumn]>dualTolerance) {
     
    25032618          smallerSlack=newSlack;
    25042619        }
    2505         double infeasibility = fabs(newPrimal+upperSlack[iColumn]-upper[iColumn]);
    2506         if (infeasibility>maximumBoundInfeasibility) {
    2507           maximumBoundInfeasibility=infeasibility;
    2508         }
     2620        if (!fakeUpper(iColumn)) {
     2621          double infeasibility = fabs(newPrimal+upperSlack[iColumn]-upper[iColumn]);
     2622          if (infeasibility>maximumBoundInfeasibility) {
     2623            maximumBoundInfeasibility=infeasibility;
     2624          }
     2625        }
    25092626        if (upperSlack[iColumn]<=kill&&fabs(newPrimal-upper[iColumn])<=kill) {
    25102627          //may be better to leave at value?
     
    25162633        }
    25172634      }
     2635      if ((!lowerBound(iColumn))&&
     2636               (!upperBound(iColumn))) {
     2637        double gap;
     2638        if (fabs(newPrimal)>eFree) {
     2639          gap=fabs(newPrimal);
     2640        } else {
     2641          gap=eFree;
     2642        }
     2643        zVec[iColumn]=gap*freeMultiplier;
     2644        tVec[iColumn]=zVec[iColumn];
     2645        //fake to give correct result
     2646        s=1.0;
     2647        t=1.0;
     2648        tValue=0.0;
     2649        // If dual infeasible then adjust more
     2650        double dualInfeasibility=reducedCost-zVec[iColumn]+tVec[iColumn];
     2651        if (fabs(dualInfeasibility)>2.0*dualTolerance&&
     2652            fabs(dualInfeasibility)>0.1*sumDualInfeasibilities_) {
     2653          zValue=2.0*freeMultiplier+max(qDiagonal,1.0e-4);
     2654          zValue=2.0*freeMultiplier+max(qDiagonal*1.0e-5,1.0e-4);
     2655          zValue=2.0*freeMultiplier+qDiagonal;
     2656        } else {
     2657          zValue=2.0*freeMultiplier+qDiagonal;
     2658        }
     2659        // take out of infeasibility
     2660        //reducedCost=zVec[iColumn]-tVec[iColumn];
     2661      }
    25182662      primal[iColumn]=newPrimal;
    25192663      if (fabs(newPrimal)>solutionNorm) {
     
    25482692        double dualInfeasibility=reducedCost-zVec[iColumn]+tVec[iColumn];
    25492693        if (fabs(dualInfeasibility)>dualTolerance) {
    2550           dualFake+=newPrimal*dualInfeasibility;
     2694          if ((!lowerBound(iColumn))&&
     2695                 (!upperBound(iColumn))) {
     2696            dualFake+=newPrimal*dualInfeasibility;
     2697          } else {
     2698            dualFake+=newPrimal*dualInfeasibility;
     2699          }
    25512700        }
    25522701        dualInfeasibility=fabs(dualInfeasibility);
     
    27532902    double w4=-deltaT[iColumn]*deltaX[iColumn];
    27542903    if (lowerBound(iColumn)) {
    2755       w3+=deltaZ[iColumn]*(primal[iColumn]-lowerSlack[iColumn]-lower[iColumn]);
     2904      if (!fakeLower(iColumn)) {
     2905        w3+=deltaZ[iColumn]*(primal[iColumn]-lowerSlack[iColumn]-lower[iColumn]);
     2906      }
    27562907      product+=w3;
    27572908    }
    27582909    if (upperBound(iColumn)) {
    2759       w4+=deltaT[iColumn]*(-primal[iColumn]-upperSlack[iColumn]+upper[iColumn]);
     2910      if (!fakeUpper(iColumn)) {
     2911        w4+=deltaT[iColumn]*(-primal[iColumn]-upperSlack[iColumn]+upper[iColumn]);
     2912      }
    27602913      product+=w4;
    27612914    }
  • trunk/ClpPresolve.cpp

    r368 r378  
    721721
    722722  if (prob.colstat_)
    723     prob.check_nbasic();
     723   prob.check_nbasic();
    724724 
    725725#if     DEBUG_PRESOLVE
  • trunk/ClpSolve.cpp

    r377 r378  
    895895  } else if (method==ClpSolve::useBarrier) {
    896896    //printf("***** experimental pretty crude barrier\n");
    897     //#define BORROW
     897    //#define SAVEIT 1
     898#ifndef SAVEIT
     899#define BORROW
     900#endif
    898901#ifdef BORROW
    899902    ClpInterior barrier;
     
    913916    barrier.setCholesky(cholesky);
    914917#endif
    915     //#define SAVEIT 1
     918    int numberRows = model2->numberRows();
     919    int numberColumns = model2->numberColumns();
    916920#ifndef SAVEIT
    917921    barrier.primalDual();
     
    924928    // move solutions
    925929    CoinMemcpyN(model2->primalRowSolution(),
    926                 barrier.numberRows(),barrier.primalRowSolution());
     930                numberRows,barrier.primalRowSolution());
    927931    CoinMemcpyN(model2->dualRowSolution(),
    928                 barrier.numberRows(),barrier.dualRowSolution());
     932                numberRows,barrier.dualRowSolution());
    929933    CoinMemcpyN(model2->primalColumnSolution(),
    930                 barrier.numberColumns(),barrier.primalColumnSolution());
     934                numberColumns,barrier.primalColumnSolution());
    931935    CoinMemcpyN(model2->dualColumnSolution(),
    932                 barrier.numberColumns(),barrier.dualColumnSolution());
     936                numberColumns,barrier.dualColumnSolution());
    933937#endif
    934938    time2 = CoinCpuTime();
     
    938942      <<CoinMessageEol;
    939943    timeX=time2;
    940     if (barrier.maximumBarrierIterations()&&barrier.status()<4) {
     944    int maxIts = barrier.maximumBarrierIterations();
     945    int barrierStatus=barrier.status();
     946    double gap = barrier.complementarityGap();
     947    // could get which variables are fixed - would help crossover *****
     948#ifdef BORROW   
     949    barrier.returnModel(*model2);
     950    double * rowPrimal = new double [numberRows];
     951    double * columnPrimal = new double [numberColumns];
     952    double * rowDual = new double [numberRows];
     953    double * columnDual = new double [numberColumns];
     954    // move solutions other way
     955    CoinMemcpyN(model2->primalRowSolution(),
     956                numberRows,rowPrimal);
     957    CoinMemcpyN(model2->dualRowSolution(),
     958                numberRows,rowDual);
     959    CoinMemcpyN(model2->primalColumnSolution(),
     960                numberColumns,columnPrimal);
     961    CoinMemcpyN(model2->dualColumnSolution(),
     962                  numberColumns,columnDual);
     963#else
     964    double * rowPrimal = barrier.primalRowSolution();
     965    double * columnPrimal = barrier.primalColumnSolution();
     966    double * rowDual = barrier.dualRowSolution();
     967    double * columnDual = barrier.dualColumnSolution();
     968    // move solutions
     969    CoinMemcpyN(rowPrimal,
     970                numberRows,model2->primalRowSolution());
     971    CoinMemcpyN(rowDual,
     972                numberRows,model2->dualRowSolution());
     973    CoinMemcpyN(columnPrimal,
     974                numberColumns,model2->primalColumnSolution());
     975    CoinMemcpyN(columnDual,
     976                  numberColumns,model2->dualColumnSolution());
     977#endif
     978    if (maxIts&&barrierStatus<4) {
    941979      printf("***** crossover - needs more thought on difficult models\n");
    942       // move solutions
    943       CoinMemcpyN(barrier.primalRowSolution(),
    944                   model2->numberRows(),model2->primalRowSolution());
    945       CoinMemcpyN(barrier.dualRowSolution(),
    946                   model2->numberRows(),model2->dualRowSolution());
    947       CoinMemcpyN(barrier.primalColumnSolution(),
    948                   model2->numberColumns(),model2->primalColumnSolution());
    949       CoinMemcpyN(barrier.dualColumnSolution(),
    950                   model2->numberColumns(),model2->dualColumnSolution());
    951980#if SAVEIT==1
    952981      model2->ClpSimplex::saveModel("xx.save");
     
    959988      // throw some into basis
    960989      {
    961         int numberColumns = model2->numberColumns();
    962         int numberRows = model2->numberRows();
    963990        double * dsort = new double[numberColumns];
    964991        int * sort = new int[numberColumns];
     
    9951022        delete [] dsort;
    9961023      }
    997       // just primal values pass
    998       model2->primal(2);
    999       // move solutions
    1000       CoinMemcpyN(model2->primalRowSolution(),
    1001                   model2->numberRows(),barrier.primalRowSolution());
    1002       CoinMemcpyN(barrier.dualRowSolution(),
    1003                   model2->numberRows(),model2->dualRowSolution());
    1004       CoinMemcpyN(model2->primalColumnSolution(),
    1005                   model2->numberColumns(),barrier.primalColumnSolution());
    1006       CoinMemcpyN(barrier.dualColumnSolution(),
    1007                   model2->numberColumns(),model2->dualColumnSolution());
    1008       //model2->primal(1);
    1009       // clean up reduced costs and flag variables
    1010       {
    1011         int numberColumns = model2->numberColumns();
    1012         double * dj = model2->dualColumnSolution();
    1013         double * cost = model2->objective();
    1014         double * saveCost = new double[numberColumns];
    1015         memcpy(saveCost,cost,numberColumns*sizeof(double));
    1016         double * saveLower = new double[numberColumns];
    1017         double * lower = model2->columnLower();
    1018         memcpy(saveLower,lower,numberColumns*sizeof(double));
    1019         double * saveUpper = new double[numberColumns];
    1020         double * upper = model2->columnUpper();
    1021         memcpy(saveUpper,upper,numberColumns*sizeof(double));
    1022         int i;
    1023         double tolerance = 10.0*dualTolerance_;
    1024         for ( i=0;i<numberColumns;i++) {
    1025           if (model2->getStatus(i)==basic) {
    1026             dj[i]=0.0;
    1027           } else if (model2->getStatus(i)==atLowerBound) {
    1028             if (optimizationDirection_*dj[i]<tolerance) {
    1029               if (optimizationDirection_*dj[i]<0.0) {
    1030                 //if (dj[i]<-1.0e-3)
    1031                 //printf("bad dj at lb %d %g\n",i,dj[i]);
    1032                 cost[i] -= dj[i];
    1033                 dj[i]=0.0;
     1024      if (gap<1.0e-3*((double) (numberRows+numberColumns))) {
     1025        // just primal values pass
     1026        model2->primal(2);
     1027        // save primal solution and copy back dual
     1028        CoinMemcpyN(model2->primalRowSolution(),
     1029                    numberRows,rowPrimal);
     1030        CoinMemcpyN(rowDual,
     1031                    numberRows,model2->dualRowSolution());
     1032        CoinMemcpyN(model2->primalColumnSolution(),
     1033                    numberColumns,columnPrimal);
     1034        CoinMemcpyN(columnDual,
     1035                    numberColumns,model2->dualColumnSolution());
     1036        //model2->primal(1);
     1037        // clean up reduced costs and flag variables
     1038        {
     1039          double * dj = model2->dualColumnSolution();
     1040          double * cost = model2->objective();
     1041          double * saveCost = new double[numberColumns];
     1042          memcpy(saveCost,cost,numberColumns*sizeof(double));
     1043          double * saveLower = new double[numberColumns];
     1044          double * lower = model2->columnLower();
     1045          memcpy(saveLower,lower,numberColumns*sizeof(double));
     1046          double * saveUpper = new double[numberColumns];
     1047          double * upper = model2->columnUpper();
     1048          memcpy(saveUpper,upper,numberColumns*sizeof(double));
     1049          int i;
     1050          double tolerance = 10.0*dualTolerance_;
     1051          for ( i=0;i<numberColumns;i++) {
     1052            if (model2->getStatus(i)==basic) {
     1053              dj[i]=0.0;
     1054            } else if (model2->getStatus(i)==atLowerBound) {
     1055              if (optimizationDirection_*dj[i]<tolerance) {
     1056                if (optimizationDirection_*dj[i]<0.0) {
     1057                  //if (dj[i]<-1.0e-3)
     1058                  //printf("bad dj at lb %d %g\n",i,dj[i]);
     1059                  cost[i] -= dj[i];
     1060                  dj[i]=0.0;
     1061                }
     1062              } else {
     1063                upper[i]=lower[i];
    10341064              }
    1035             } else {
    1036               upper[i]=lower[i];
    1037             }
    1038           } else if (model2->getStatus(i)==atUpperBound) {
    1039             if (optimizationDirection_*dj[i]>tolerance) {
    1040               if (optimizationDirection_*dj[i]>0.0) {
    1041                 //if (dj[i]>1.0e-3)
    1042                 //printf("bad dj at ub %d %g\n",i,dj[i]);
    1043                 cost[i] -= dj[i];
    1044                 dj[i]=0.0;
     1065            } else if (model2->getStatus(i)==atUpperBound) {
     1066              if (optimizationDirection_*dj[i]>tolerance) {
     1067                if (optimizationDirection_*dj[i]>0.0) {
     1068                  //if (dj[i]>1.0e-3)
     1069                  //printf("bad dj at ub %d %g\n",i,dj[i]);
     1070                  cost[i] -= dj[i];
     1071                  dj[i]=0.0;
     1072                }
     1073              } else {
     1074                lower[i]=upper[i];
    10451075              }
    1046             } else {
    1047               lower[i]=upper[i];
    10481076            }
    10491077          }
    1050         }
    1051         // just dual values pass
    1052         //model2->setLogLevel(63);
    1053         //model2->setFactorizationFrequency(1);
    1054         model2->dual(2);
    1055         memcpy(cost,saveCost,numberColumns*sizeof(double));
    1056         delete [] saveCost;
    1057         memcpy(lower,saveLower,numberColumns*sizeof(double));
    1058         delete [] saveLower;
    1059         memcpy(upper,saveUpper,numberColumns*sizeof(double));
    1060         delete [] saveUpper;
    1061       }
    1062       // and finish
    1063       // move solutions
    1064       CoinMemcpyN(barrier.primalRowSolution(),
    1065                   model2->numberRows(),model2->primalRowSolution());
    1066       CoinMemcpyN(barrier.primalColumnSolution(),
    1067                   model2->numberColumns(),model2->primalColumnSolution());
     1078          // just dual values pass
     1079          //model2->setLogLevel(63);
     1080          //model2->setFactorizationFrequency(1);
     1081          model2->dual(2);
     1082          memcpy(cost,saveCost,numberColumns*sizeof(double));
     1083          delete [] saveCost;
     1084          memcpy(lower,saveLower,numberColumns*sizeof(double));
     1085          delete [] saveLower;
     1086          memcpy(upper,saveUpper,numberColumns*sizeof(double));
     1087          delete [] saveUpper;
     1088        }
     1089        // and finish
     1090        // move solutions
     1091        CoinMemcpyN(rowPrimal,
     1092                    numberRows,model2->primalRowSolution());
     1093        CoinMemcpyN(columnPrimal,
     1094                    numberColumns,model2->primalColumnSolution());
     1095      }
    10681096      model2->primal(1);
    10691097#else
     
    10711099      model2->primal(1);
    10721100#endif
    1073     } else if (barrier.status()==4) {
     1101    } else if (barrierStatus==4) {
    10741102      // memory problems
    10751103      model2->setPerturbation(savePerturbation);
     
    10781106    }
    10791107#ifdef BORROW
    1080     barrier.returnModel(*model2);
     1108    delete [] rowPrimal;
     1109    delete [] columnPrimal;
     1110    delete [] rowDual;
     1111    delete [] columnDual;
    10811112#endif
    10821113    model2->setPerturbation(savePerturbation);
  • trunk/include/ClpInterior.hpp

    r369 r378  
    137137              bool keepNames=false,
    138138              bool ignoreErrors = false);
     139  /** Borrow model.  This is so we dont have to copy large amounts
     140      of data around.  It assumes a derived class wants to overwrite
     141      an empty model with a real one - while it does an algorithm.
     142      This is same as ClpModel one. */
     143  void borrowModel(ClpModel & otherModel);
     144  /** Return model - updates any scalars */
     145  void returnModel(ClpModel & otherModel);
    139146  //@}
    140147
     
    182189  inline void setDiagonalPerturbation(double value)
    183190  { diagonalPerturbation_=value;};
     191  /// ComplementarityGap
     192  inline double complementarityGap() const
     193          { return complementarityGap_;} ;
    184194  //@}
    185195
Note: See TracChangeset for help on using the changeset viewer.