Changeset 1676


Ignore:
Timestamp:
Jan 20, 2011 12:00:28 PM (9 years ago)
Author:
forrest
Message:

For GUB

Location:
trunk/Clp/src
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpDynamicExampleMatrix.cpp

    r1665 r1676  
    655655          }
    656656     }
    657      ClpDynamicMatrix::createVariable(model, bestSequence);
     657     ClpDynamicMatrix::createVariable(model, bestSequence/*, bestSequence2*/);
    658658     // clear for next iteration
    659659     savedBestSequence_ = -1;
  • trunk/Clp/src/ClpDynamicMatrix.cpp

    r1665 r1676  
    9494     lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
    9595     upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
    96      status_ = ClpCopyOfArray(rhs.status_, numberSets_);
     96     status_ = ClpCopyOfArray(rhs.status_, 2*numberSets_);
    9797     model_ = rhs.model_;
    9898     sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
     
    110110     maximumGubColumns_ = rhs.maximumGubColumns_;
    111111     maximumElements_ = rhs.maximumElements_;
    112      startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_);
     112     startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_+1);
    113113     next_ = ClpCopyOfArray(rhs.next_, maximumGubColumns_);
    114114     startColumn_ = ClpCopyOfArray(rhs.startColumn_, maximumGubColumns_ + 1);
     
    119119     columnLower_ = ClpCopyOfArray(rhs.columnLower_, maximumGubColumns_);
    120120     columnUpper_ = ClpCopyOfArray(rhs.columnUpper_, maximumGubColumns_);
    121      dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, maximumGubColumns_);
     121     dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, 2*maximumGubColumns_);
    122122}
    123123
     
    143143     else
    144144          maximumElements_ = 0;
    145      startSet_ = new int [numberSets_];
     145     startSet_ = new int [numberSets_+1];
    146146     next_ = new int [maximumGubColumns_];
    147147     // fill in startSet and next
     
    157157          }
    158158     }
     159     startSet_[numberSets_] = starts[numberSets_];
    159160     int numberColumns = model->numberColumns();
    160161     int numberRows = model->numberRows();
     
    165166     int frequency = model->factorizationFrequency();
    166167     int numberGubInSmall = numberRows + frequency + CoinMin(frequency, numberSets_) + 4;
     168     // But we may have two per row + one for incoming (make it two)
     169     numberGubInSmall = CoinMax(2*numberRows+2,numberGubInSmall);
    167170     // for small problems this could be too big
    168171     //numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_);
     
    256259     keyVariable_ = new int[numberSets_];
    257260     if (status) {
    258           status_ = ClpCopyOfArray(status, numberSets_);
     261          status_ = ClpCopyOfArray(status, 2*numberSets_);
    259262          assert (dynamicStatus);
    260           dynamicStatus_ = ClpCopyOfArray(dynamicStatus, numberGubColumns_);
     263          dynamicStatus_ = ClpCopyOfArray(dynamicStatus, 2*numberGubColumns_);
    261264     } else {
    262           status_ = new unsigned char [numberSets_];
     265          status_ = new unsigned char [2*numberSets_];
    263266          memset(status_, 0, numberSets_);
    264267          int i;
     
    267270               setStatus(i, ClpSimplex::basic);
    268271          }
    269           dynamicStatus_ = new unsigned char [numberGubColumns_];
     272          dynamicStatus_ = new unsigned char [2*numberGubColumns_];
    270273          memset(dynamicStatus_, 0, numberGubColumns_); // for clarity
    271274          for (i = 0; i < numberGubColumns_; i++)
     
    353356          lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
    354357          upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
    355           status_ = ClpCopyOfArray(rhs.status_, numberSets_);
     358          status_ = ClpCopyOfArray(rhs.status_, 2*numberSets_);
    356359          model_ = rhs.model_;
    357360          sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
     
    369372          maximumGubColumns_ = rhs.maximumGubColumns_;
    370373          maximumElements_ = rhs.maximumElements_;
    371           startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_);
     374          startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_+1);
    372375          next_ = ClpCopyOfArray(rhs.next_, maximumGubColumns_);
    373376          startColumn_ = ClpCopyOfArray(rhs.startColumn_, maximumGubColumns_ + 1);
     
    378381          columnLower_ = ClpCopyOfArray(rhs.columnLower_, maximumGubColumns_);
    379382          columnUpper_ = ClpCopyOfArray(rhs.columnUpper_, maximumGubColumns_);
    380           dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, maximumGubColumns_);
     383          dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, 2*maximumGubColumns_);
    381384     }
    382385     return *this;
     
    533536                                   // check flagged variable and correct dj
    534537                                   if (!flagged(iSequence)) {
     538                                     if (false/*status == atLowerBound
     539                                                &&keyValue(iSet)<1.0e-7*/) {
     540                                       // can't come in because
     541                                       // of ones at ub
     542                                       numberWanted++;
     543                                     } else {
     544                                     
    535545                                        bestDj = value;
    536546                                        bestSequence = structuralOffset + iSequence;
    537547                                        bestDjMod = djMod;
    538548                                        bestSet = iSet;
     549                                     }
    539550                                   } else {
    540551                                        // just to make sure we don't exit before got something
     
    576587                           )
    577588{
    578      //forceRefresh=true;
     589     // forceRefresh=true;printf("take out forceRefresh\n");
    579590     if (!model_->numberIterations())
    580591          forceRefresh = true;
     
    735746                                   } else if (getDynamicStatus(j) == atUpperBound) {
    736747                                        value = columnUpper_[j];
     748                                        assert (value<1.0e30);
    737749                                   } else if (getDynamicStatus(j) == soloKey) {
    738750                                        value = keyValue(iSet);
    739                                    }
     751                                   }
    740752                                   objectiveOffset += value * cost_[j];
    741753                              }
     
    753765                    for (iSet = 0; iSet < numberSets_; iSet++) {
    754766                         int kRow = toIndex_[iSet];
    755                          if (kRow >= 0)
    756                               kRow += numberStaticRows_;
    757                          int j = startSet_[iSet];
    758                          while (j >= 0) {
    759                               double value = solution[j];
    760                               if (value) {
    761                                    for (CoinBigIndex k = startColumn_[j]; k < startColumn_[j+1]; k++) {
    762                                         int iRow = row_[k];
    763                                         rhsOffset_[iRow] -= element_[k] * value;
    764                                    }
    765                                    if (kRow >= 0)
    766                                         rhsOffset_[kRow] -= value;
    767                               }
    768                               j = next_[j]; //onto next in set
    769                          }
     767                         if (kRow >= 0)
     768                           kRow += numberStaticRows_;
     769                         int j = startSet_[iSet];
     770                         while (j >= 0) {
     771                           //? use DynamicStatus status = getDynamicStatus(j);
     772                           double value = solution[j];
     773                           if (value) {
     774                             for (CoinBigIndex k = startColumn_[j]; k < startColumn_[j+1]; k++) {
     775                               int iRow = row_[k];
     776                               rhsOffset_[iRow] -= element_[k] * value;
     777                             }
     778                             if (kRow >= 0)
     779                               rhsOffset_[kRow] -= value;
     780                           }
     781                           j = next_[j]; //onto next in set
     782                         }
    770783                    }
    771784                    delete [] solution;
     
    10731086     // save status
    10741087     case 5: {
     1088       memcpy(status_+numberSets_,status_,numberSets_);
     1089       memcpy(dynamicStatus_+maximumGubColumns_,
     1090              dynamicStatus_,maximumGubColumns_);
    10751091     }
    10761092     break;
    10771093     // restore status
    10781094     case 6: {
    1079           // maybe mismatch - redo all in small model
    1080           //for (int i=firstDynamic_;i<firstAvailable_;i++) {
    1081           //int jColumn = id_[i-firstDynamic_];
    1082           //setDynamicStatus(jColumn,inSmall);
    1083           //}
     1095       memcpy(status_,status_+numberSets_,numberSets_);
     1096       memcpy(dynamicStatus_,dynamicStatus_+maximumGubColumns_,
     1097              maximumGubColumns_);
     1098       initialProblem();
    10841099     }
    10851100     break;
     
    11721187          // first flag
    11731188          if (number >= firstDynamic_ && number < lastDynamic_) {
    1174                assert (number == model->sequenceIn());
    1175                // id will be sitting at firstAvailable
    1176                int sequence = id_[firstAvailable_-firstDynamic_];
    1177                setFlagged(sequence);
    1178                model->clearFlagged(firstAvailable_);
    1179           }
     1189            int sequence = id_[number-firstDynamic_];
     1190            setFlagged(sequence);
     1191            //model->clearFlagged(number);
     1192          } else if (number>=model_->numberColumns()+numberStaticRows_) {
     1193            // slack
     1194            int iSet = fromIndex_[number-model_->numberColumns()-
     1195                                  numberStaticRows_];
     1196            setFlaggedSlack(iSet);
     1197            //model->clearFlagged(number);
     1198          }
    11801199     }
    11811200     case 11: {
    1182           int sequenceIn = model->sequenceIn();
    1183           if (sequenceIn >= firstDynamic_ && sequenceIn < lastDynamic_) {
    1184                assert (number == model->sequenceIn());
     1201          //int sequenceIn = model->sequenceIn();
     1202          if (number >= firstDynamic_ && number < lastDynamic_) {
     1203            //assert (number == model->sequenceIn());
    11851204               // take out variable (but leave key)
    11861205               double * cost = model->costRegion();
     
    11901209               int * length = matrix_->getMutableVectorLengths();
    11911210#ifndef NDEBUG
    1192                if (length[sequenceIn]) {
     1211               if (length[number]) {
    11931212                    int * row = matrix_->getMutableIndices();
    11941213                    CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
    1195                     int iRow = row[startColumn[sequenceIn] + length[sequenceIn] - 1];
     1214                    int iRow = row[startColumn[number] + length[number] - 1];
    11961215                    int which = iRow - numberStaticRows_;
    11971216                    assert (which >= 0);
     
    12101229
    12111230               // not really in small problem
    1212                int iBig = id_[sequenceIn-firstDynamic_];
    1213                if (model->getStatus(sequenceIn) == ClpSimplex::atLowerBound) {
     1231               int iBig = id_[number-firstDynamic_];
     1232               if (model->getStatus(number) == ClpSimplex::atLowerBound) {
    12141233                    setDynamicStatus(iBig, atLowerBound);
    12151234                    if (columnLower_)
    1216                          modifyOffset(sequenceIn, columnLower_[iBig]);
     1235                         modifyOffset(number, columnLower_[iBig]);
    12171236               } else {
    12181237                    setDynamicStatus(iBig, atUpperBound);
    1219                     modifyOffset(sequenceIn, columnUpper_[iBig]);
    1220                }
     1238                    modifyOffset(number, columnUpper_[iBig]);
     1239               }
     1240          } else if (number>=model_->numberColumns()+numberStaticRows_) {
     1241            // slack
     1242            int iSet = fromIndex_[number-model_->numberColumns()-
     1243                                  numberStaticRows_];
     1244            printf("what now - set %d\n",iSet);
    12211245          }
    12221246     }
     
    13021326               double value = solution[iColumn];
    13031327               bool toLowerBound = true;
     1328               assert (jColumn>=0);assert (iColumn>=0);
    13041329               if (columnUpper_) {
    13051330                    if (!columnLower_) {
    13061331                         if (fabs(value - columnUpper_[jColumn]) < fabs(value))
    13071332                              toLowerBound = false;
    1308                     } else if (fabs(value - columnUpper_[jColumn]) < fabs(value - columnLower_[iColumn])) {
     1333                    } else if (fabs(value - columnUpper_[jColumn]) < fabs(value - columnLower_[jColumn])) {
    13091334                         toLowerBound = false;
    13101335                    }
     
    15091534               n++;
    15101535          }
     1536          int k=0;
     1537          for (int j=startSet_[i];j<startSet_[i+1];j++) {
     1538            if (getDynamicStatus(j)==soloKey)
     1539              k++;
     1540          }
     1541          assert (k<2);
    15111542     }
    15121543     assert (n == numberSets_);
     
    15971628                    j = next_[j]; //onto next in set
    15981629               }
    1599           }
     1630#if 0
     1631               // slack must be feasible
     1632               double oldValue=value;
     1633               value = CoinMax(value,lowerSet_[iSet]);
     1634               value = CoinMin(value,upperSet_[iSet]);
     1635               if (value!=oldValue)
     1636                 printf("using %g (not %g) for slack on set %d (%g,%g)\n",
     1637                        value,oldValue,iSet,lowerSet_[iSet],upperSet_[iSet]);
     1638#endif
     1639            }
    16001640     }
    16011641     return value;
     
    16761716                    model->setStatus(iSequence, getStatus(savedBestSet_));
    16771717                    model->djRegion()[iSequence] = savedBestGubDual_;
    1678                     if (getStatus(savedBestSet_) == ClpSimplex::atUpperBound)
    1679                          solution[iSequence] = columnUpper[iSequence];
    1680                     else
    1681                          solution[iSequence] = columnLower[iSequence];
    1682                     // create variable and pivot in
     1718                    solution[iSequence] = valueOfKey;
     1719                    // create variable and pivot in
    16831720                    int key = keyVariable_[savedBestSet_];
    16841721                    setDynamicStatus(key, inSmall);
     
    17701807                    model->setStatus(iSequence, ClpSimplex::basic);
    17711808                    model->djRegion()[iSequence] = 0.0;
    1772                     solution[iSequence] = shift;
     1809                    solution[iSequence] = valueOfKey+shift;
    17731810                    rhsOffset_[newRow] = -shift; // sign?
    17741811               }
     
    18161853               //printf("best %d\n",bestSequence2);
    18171854               model->solutionRegion()[firstAvailable_] = 0.0;
     1855               model->clearFlagged(firstAvailable_);
    18181856               if (!columnLower_ && !columnUpper_) {
    18191857                    model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
     
    21092147          int numberActive = 0;
    21102148          int whichKey = -1;
    2111           if (getStatus(iSet) == ClpSimplex::basic)
     2149          if (getStatus(iSet) == ClpSimplex::basic) {
    21122150               whichKey = maximumGubColumns_ + iSet;
    2113           else
     2151               numberActive = 1;
     2152          } else {
    21142153               whichKey = -1;
     2154          }
    21152155          int j = startSet_[iSet];
    21162156          while (j >= 0) {
    21172157               assert (getDynamicStatus(j) != soloKey || whichKey == -1);
    2118                if (getDynamicStatus(j) == inSmall)
     2158               if (getDynamicStatus(j) == inSmall) {
    21192159                    numberActive++;
    2120                else if (getDynamicStatus(j) == soloKey)
     2160               } else if (getDynamicStatus(j) == soloKey) {
    21212161                    whichKey = j;
     2162                    numberActive++;
     2163               }
    21222164               j = next_[j]; //onto next in set
    21232165          }
    2124           if (getStatus(iSet) == ClpSimplex::basic && numberActive)
    2125                numberActive++;
    2126           if (numberActive) {
    2127                assert (numberActive > 1);
     2166          if (numberActive > 1) {
    21282167               int iRow = numberActiveSets_ + numberStaticRows_;
    21292168               rowSolution[iRow] = 0.0;
     
    21892228                              else
    21902229                                   columnUpper[firstAvailable_] = COIN_DBL_MAX;
    2191                               if (status == atLowerBound) {
     2230                              if (status != atUpperBound) {
    21922231                                   solution[firstAvailable_] = columnLower[firstAvailable_];
    21932232                              } else {
     
    22032242               toIndex_[iSet] = numberActiveSets_;
    22042243               fromIndex_[numberActiveSets_++] = iSet;
     2244          } else {
     2245            // solo key
     2246            bool needKey;
     2247            if (numberActive) {
     2248              if (whichKey<maximumGubColumns_) {
     2249                // structural - assume ok
     2250                needKey = false;
     2251              } else {
     2252                // slack
     2253                keyVariable_[iSet] = maximumGubColumns_ + iSet;
     2254                double value = keyValue(iSet);
     2255                if (value<lowerSet_[iSet]-1.0e-8||
     2256                    value>upperSet_[iSet]+1.0e-8)
     2257                  needKey=true;
     2258              }
     2259            } else {
     2260              needKey = true;
     2261            }
     2262            if (needKey) {
     2263              // all to lb then up some (slack/null if possible)
     2264              int length=99999999;
     2265              int which=-1;
     2266              double sum=0.0;
     2267              for (int iColumn=startSet_[iSet];iColumn<startSet_[iSet+1];iColumn++) {
     2268                setDynamicStatus(iColumn,atLowerBound);
     2269                sum += columnLower_[iColumn];
     2270                if (length>startColumn_[iColumn+1]-startColumn_[iColumn]) {
     2271                  which=iColumn;
     2272                  length=startColumn_[iColumn+1]-startColumn_[iColumn];
     2273                }
     2274              }
     2275              if (sum>lowerSet_[iSet]-1.0e-8) {
     2276                // slack can be basic
     2277                setStatus(iSet,ClpSimplex::basic);
     2278                keyVariable_[iSet] = maximumGubColumns_ + iSet;
     2279              } else {
     2280                // use shortest
     2281                setDynamicStatus(which,soloKey);
     2282                keyVariable_[iSet] = which;
     2283                setStatus(iSet,ClpSimplex::atLowerBound);
     2284              }
     2285            }
    22052286          }
    22062287          assert (toIndex_[iSet] >= 0 || whichKey >= 0);
    22072288          keyVariable_[iSet] = whichKey;
    22082289     }
     2290     numberActiveColumns_ = firstAvailable_;
    22092291     return;
     2292}
     2293// Writes out model (without names)
     2294void
     2295ClpDynamicMatrix::writeMps(const char * name)
     2296{
     2297  int numberTotalRows = numberStaticRows_+numberSets_;
     2298  int numberTotalColumns = firstDynamic_+numberGubColumns_;
     2299  // over estimate
     2300  int numberElements = getNumElements()+startColumn_[numberGubColumns_]
     2301    + numberGubColumns_;
     2302  double * columnLower = new double [numberTotalColumns];
     2303  double * columnUpper = new double [numberTotalColumns];
     2304  double * cost = new double [numberTotalColumns];
     2305  double * rowLower = new double [numberTotalRows];
     2306  double * rowUpper = new double [numberTotalRows];
     2307  CoinBigIndex * start = new CoinBigIndex[numberTotalColumns+1];
     2308  int * row = new int [numberElements];
     2309  double * element = new double [numberElements];
     2310  // Fill in
     2311  const CoinBigIndex * startA = getVectorStarts();
     2312  const int * lengthA = getVectorLengths();
     2313  const int * rowA = getIndices();
     2314  const double * elementA = getElements();
     2315  const double * columnLowerA = model_->columnLower();
     2316  const double * columnUpperA = model_->columnUpper();
     2317  const double * costA = model_->objective();
     2318  const double * rowLowerA = model_->rowLower();
     2319  const double * rowUpperA = model_->rowUpper();
     2320  start[0]=0;
     2321  numberElements=0;
     2322  for (int i=0;i<firstDynamic_;i++) {
     2323    columnLower[i] = columnLowerA[i];
     2324    columnUpper[i] = columnUpperA[i];
     2325    cost[i] = costA[i];
     2326    for (CoinBigIndex j = startA[i];j<startA[i]+lengthA[i];j++) {
     2327      row[numberElements] = rowA[j];
     2328      element[numberElements++]=elementA[j];
     2329    }
     2330    start[i+1]=numberElements;
     2331  }
     2332  for (int i=0;i<numberStaticRows_;i++) {
     2333    rowLower[i] = rowLowerA[i];
     2334    rowUpper[i] = rowUpperA[i];
     2335  }
     2336  int putC=firstDynamic_;
     2337  int putR=numberStaticRows_;
     2338  for (int i=0;i<numberSets_;i++) {
     2339    rowLower[putR]=lowerSet_[i];
     2340    rowUpper[putR]=upperSet_[i];
     2341    for (CoinBigIndex k=startSet_[i];k<startSet_[i+1];k++) {
     2342      columnLower[putC]=columnLower_[k];
     2343      columnUpper[putC]=columnUpper_[k];
     2344      cost[putC]=cost_[k];
     2345      putC++;
     2346      for (CoinBigIndex j = startColumn_[k];j<startColumn_[k+1];j++) {
     2347        row[numberElements] = row_[j];
     2348        element[numberElements++]=element_[j];
     2349      }
     2350      row[numberElements] = putR;
     2351      element[numberElements++]=1.0;
     2352      start[putC]=numberElements;
     2353    }
     2354    putR++;
     2355  }
     2356
     2357  assert (putR==numberTotalRows);
     2358  assert (putC==numberTotalColumns);
     2359  ClpSimplex modelOut;
     2360  modelOut.loadProblem(numberTotalColumns,numberTotalRows,
     2361                       start,row,element,
     2362                       columnLower,columnUpper,cost,
     2363                       rowLower,rowUpper);
     2364  modelOut.writeMps(name);
     2365  delete [] columnLower;
     2366  delete [] columnUpper;
     2367  delete [] cost;
     2368  delete [] rowLower;
     2369  delete [] rowUpper;
     2370  delete [] start;
     2371  delete [] row;
     2372  delete [] element;
    22102373}
    22112374// Adds in a column to gub structure (called from descendant)
     
    22362399                    if (columnUpper_ && upper != columnUpper_[j])
    22372400                         odd = true;
    2238                     if (odd)
     2401                    if (odd) {
    22392402                         printf("seems odd - same els but cost,lo,up are %g,%g,%g and %g,%g,%g\n",
    22402403                                cost, lower, upper, cost_[j],
    22412404                                columnLower_ ? columnLower_[j] : 0.0,
    22422405                                columnUpper_ ? columnUpper_[j] : 1.0e100);
    2243                     else
     2406                    } else {
     2407                         setDynamicStatus(j, status);
    22442408                         return j;
     2409                    }
    22452410               }
    22462411          }
     
    22592424          for (i = 0; i < numberGubColumns_; i++) {
    22602425               CoinBigIndex end = startColumn_[i+1];
     2426               // what about ubs if column generation?
    22612427               if (getDynamicStatus(i) != atLowerBound) {
    22622428                    // keep in
  • trunk/Clp/src/ClpDynamicMatrix.hpp

    r1665 r1676  
    9595     /// Does gub crash
    9696     void gubCrash();
     97     /// Writes out model (without names)
     98     void writeMps(const char * name);
    9799     /// Populates initial matrix from dynamic status
    98100     void initialProblem();
     
    166168          st_byte = static_cast<unsigned char>(st_byte | status);
    167169     }
     170     /// Whether flagged slack
     171     inline bool flaggedSlack(int i) const {
     172          return (status_[i] & 8) != 0;
     173     }
     174     inline void setFlaggedSlack(int i) {
     175          status_[i] = static_cast<unsigned char>(status_[i] | 8);
     176     }
     177     inline void unsetFlaggedSlack(int i) {
     178          status_[i] = static_cast<unsigned char>(status_[i] & ~8);
     179     }
    168180     /// Number of sets (dynamic rows)
    169181     inline int numberSets() const {
    170182          return numberSets_;
    171183     }
     184     /// Number of possible gub variables
     185     inline int numberGubEntries() const
     186     { return startSet_[numberSets_];}
     187     /// Sets
     188     inline int * startSets() const
     189     { return startSet_;}
    172190     /// Whether flagged
    173191     inline bool flagged(int i) const {
  • trunk/Clp/src/ClpNonLinearCost.cpp

    r1665 r1676  
    608608                    }
    609609               }
     610               //#define PRINT_DETAIL7 2
     611#if PRINT_DETAIL7>1
     612               printf("NL %d sol %g bounds %g %g\n",
     613                      iSequence,solution[iSequence],
     614                      lowerValue,upperValue);
     615#endif
    610616               switch(status) {
    611617
     
    625631                                        printf("nonlincostb %d %g %g %g\n",
    626632                                               iSequence, lowerValue, solution[iSequence], lower_[iRange+2]);
     633#endif
     634#if PRINT_DETAIL7
     635                                   printf("**NL %d sol %g below %g\n",
     636                                          iSequence,solution[iSequence],
     637                                          lowerValue);
    627638#endif
    628639                                   sumInfeasibilities_ += value;
     
    642653                                        printf("nonlincostu %d %g %g %g\n",
    643654                                               iSequence, lower_[iRange-1], solution[iSequence], upperValue);
     655#endif
     656#if PRINT_DETAIL7
     657                                   printf("**NL %d sol %g above %g\n",
     658                                          iSequence,solution[iSequence],
     659                                          upperValue);
    644660#endif
    645661                                   sumInfeasibilities_ += value;
  • trunk/Clp/src/ClpSimplex.cpp

    r1665 r1676  
    19331933     int cycle = progress_.cycle(in, out,
    19341934                                 directionIn_, directionOut_);
    1935      if (cycle > 0 && objective_->type() < 2) {
     1935     if (cycle > 0 && objective_->type() < 2 && matrix_->type() < 15) {
    19361936          //if (cycle>0) {
    19371937          if (handler_->logLevel() >= 63)
     
    19451945               forceFactorization_ = CoinMax(1, cycle - off[extra]);
    19461946          } else {
    1947                // need to reject something
     1947            /* need to reject something
     1948               should be better if don't reject incoming
     1949               as it is in basis */
    19481950               int iSequence;
    1949                if (algorithm_ > 0)
    1950                     iSequence = sequenceIn_;
    1951                else
     1951               //if (algorithm_ > 0)
     1952               //   iSequence = sequenceIn_;
     1953               //else
    19521954                    iSequence = sequenceOut_;
    19531955               char x = isColumn(iSequence) ? 'C' : 'R';
     
    19931995               forceFactorization_ = -1; //off
    19941996          return 1;
    1995      } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2))) {
     1997     } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && matrix_->type()<15) {
    19961998          double random = randomNumberGenerator_.randomDouble();
    19971999          int maxNumber = (forceFactorization_ < 0) ? maximumPivots : CoinMin(forceFactorization_, maximumPivots);
     
    82508252          whatsChanged_ &= ~0xffff;
    82518253     }
     8254     double saveObjValue = objectiveValue_;
    82528255     deleteRim(getRidOfData);
     8256     if (matrix_->type()>=15)
     8257       objectiveValue_ = saveObjValue;
    82538258     // Skip message if changing algorithms
    82548259     if (problemStatus_ != 10) {
  • trunk/Clp/src/ClpSimplexOther.cpp

    r1665 r1676  
    1616#include "ClpFactorization.hpp"
    1717#include "ClpDualRowDantzig.hpp"
     18#include "ClpDynamicMatrix.hpp"
    1819#include "CoinPackedMatrix.hpp"
    1920#include "CoinIndexedVector.hpp"
     
    40684069     checkSolutionInternal();
    40694070}
     4071// Returns gub version of model or NULL
     4072ClpSimplex *
     4073ClpSimplexOther::gubVersion(int * whichRows, int * whichColumns,
     4074                            int neededGub,
     4075                            int factorizationFrequency)
     4076{
     4077  // find gub
     4078  int numberRows = this->numberRows();
     4079  int numberColumns = this->numberColumns();
     4080  int iRow, iColumn;
     4081  int * columnIsGub = new int [numberColumns];
     4082  const double * columnLower = this->columnLower();
     4083  const double * columnUpper = this->columnUpper();
     4084  int numberFixed=0;
     4085  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     4086    if (columnUpper[iColumn] == columnLower[iColumn]) {
     4087      columnIsGub[iColumn]=-2;
     4088      numberFixed++;
     4089    } else if (columnLower[iColumn]>=0) {
     4090      columnIsGub[iColumn]=-1;
     4091    } else {
     4092      columnIsGub[iColumn]=-3;
     4093    }
     4094  }
     4095  CoinPackedMatrix * matrix = this->matrix();
     4096  // get row copy
     4097  CoinPackedMatrix rowCopy = *matrix;
     4098  rowCopy.reverseOrdering();
     4099  const int * column = rowCopy.getIndices();
     4100  const int * rowLength = rowCopy.getVectorLengths();
     4101  const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
     4102  const double * element = rowCopy.getElements();
     4103  int numberNonGub = 0;
     4104  int numberEmpty = numberRows;
     4105  int * rowIsGub = new int [numberRows];
     4106  int smallestGubRow=-1;
     4107  int count=numberColumns+1;
     4108  double * rowLower = this->rowLower();
     4109  double * rowUpper = this->rowUpper();
     4110  // make sure we can get rid of upper bounds
     4111  double * fixedRow = new double [numberRows];
     4112  for (iRow = 0 ; iRow < numberRows ; iRow++) {
     4113    double sumFixed=0.0;
     4114    for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
     4115      int iColumn = column[j];
     4116      double value = columnLower[iColumn];
     4117      if (value)
     4118        sumFixed += element[j] * value;
     4119    }
     4120    fixedRow[iRow]=rowUpper[iRow]-sumFixed;
     4121  }
     4122  for (iRow = numberRows - 1; iRow >= 0; iRow--) {
     4123    bool gubRow = true;
     4124    int numberInRow=0;
     4125    double sumFixed=0.0;
     4126    double gap = fixedRow[iRow]-1.0e-12;
     4127    for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
     4128      int iColumn = column[j];
     4129      if (columnIsGub[iColumn]!=-2) {
     4130        if (element[j] != 1.0||columnIsGub[iColumn]==-3||
     4131            columnUpper[iColumn]-columnLower[iColumn]<gap) {
     4132          gubRow = false;
     4133          break;
     4134        } else {
     4135          numberInRow++;
     4136          if (columnIsGub[iColumn] >= 0) {
     4137            gubRow = false;
     4138            break;
     4139          }
     4140        }
     4141      } else {
     4142        sumFixed += columnLower[iColumn]*element[j];
     4143      }
     4144    }
     4145    if (!gubRow) {
     4146      whichRows[numberNonGub++] = iRow;
     4147      rowIsGub[iRow] = -1;
     4148    } else if (numberInRow) {
     4149      if (numberInRow<count) {
     4150        count = numberInRow;
     4151        smallestGubRow=iRow;
     4152      }
     4153      for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
     4154        int iColumn = column[j];
     4155        if (columnIsGub[iColumn]!=-2)
     4156          columnIsGub[iColumn] = iRow;
     4157      }
     4158      rowIsGub[iRow] = 0;
     4159    } else {
     4160      // empty row!
     4161      whichRows[--numberEmpty] = iRow;
     4162      rowIsGub[iRow] = -2;
     4163      if (sumFixed>rowUpper[iRow]+1.0e-4||
     4164          sumFixed<rowLower[iRow]-1.0e-4) {
     4165        fprintf(stderr,"******** No infeasible empty rows - please!\n");
     4166        abort();
     4167      }
     4168    }
     4169  }
     4170  delete [] fixedRow;
     4171  char message[100];
     4172  int numberGub = numberEmpty - numberNonGub;
     4173  if (numberGub >= neededGub) {
     4174    sprintf(message,"%d gub rows", numberGub);
     4175    handler_->message(CLP_GENERAL, messages_)
     4176      << message << CoinMessageEol;
     4177    int numberNormal = 0;
     4178    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     4179      if (columnIsGub[iColumn] < 0  && columnIsGub[iColumn] !=-2) {
     4180        whichColumns[numberNormal++] = iColumn;
     4181      }
     4182    }
     4183    if (!numberNormal) {
     4184      sprintf(message,"Putting back one gub row to make non-empty");
     4185      handler_->message(CLP_GENERAL, messages_)
     4186        << message << CoinMessageEol;
     4187      rowIsGub[smallestGubRow]=-1;
     4188      whichRows[numberNonGub++] = smallestGubRow;
     4189      for (int j = rowStart[smallestGubRow];
     4190           j < rowStart[smallestGubRow] + rowLength[smallestGubRow]; j++) {
     4191        int iColumn = column[j];
     4192        if (columnIsGub[iColumn]>=0) {
     4193          columnIsGub[iColumn]=-4;
     4194          whichColumns[numberNormal++] = iColumn;
     4195        }
     4196      }
     4197    }
     4198    std::sort(whichRows,whichRows+numberNonGub);
     4199    std::sort(whichColumns,whichColumns+numberNormal);
     4200    double * lower = CoinCopyOfArray(this->rowLower(),numberRows);
     4201    double * upper = CoinCopyOfArray(this->rowUpper(),numberRows);
     4202    // leave empty rows at end
     4203    numberEmpty = numberRows-numberEmpty;
     4204    const int * row = matrix->getIndices();
     4205    const int * columnLength = matrix->getVectorLengths();
     4206    const CoinBigIndex * columnStart = matrix->getVectorStarts();
     4207    const double * elementByColumn = matrix->getElements();
     4208    // Fixed at end
     4209    int put2 = numberColumns-numberFixed;
     4210    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     4211      if (columnIsGub[iColumn] ==-2) {
     4212        whichColumns[put2++] = iColumn;
     4213        double value = columnLower[iColumn];
     4214        for (int j = columnStart[iColumn];
     4215             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     4216          int iRow = row[j];
     4217          if (lower[iRow]>-1.0e20)
     4218            lower[iRow] -= value*element[j];
     4219          if (upper[iRow]<1.0e20)
     4220            upper[iRow] -= value*element[j];
     4221        }
     4222      }
     4223    }
     4224    int put = numberNormal;
     4225    ClpSimplex * model2 =
     4226      new ClpSimplex(this, numberNonGub, whichRows , numberNormal, whichColumns);
     4227    // scale
     4228    double * scaleArray = new double [numberRows];
     4229    for (int i=0;i<numberRows;i++) {
     4230      scaleArray[i]=1.0;
     4231      if (rowIsGub[i]==-1) {
     4232        double largest = 1.0e-30;
     4233        double smallest = 1.0e30;
     4234        for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
     4235          int iColumn = column[j];
     4236          if (columnIsGub[iColumn]!=-2) {
     4237            double value =fabs(element[j]);
     4238            largest = CoinMax(value,largest);
     4239            smallest = CoinMin(value,smallest);
     4240          }
     4241        }
     4242        double scale = CoinMax(0.001,1.0/sqrt(largest*smallest));
     4243        scaleArray[i]=scale;
     4244        if (lower[i]>-1.0e30)
     4245          lower[i] *= scale;
     4246        if (upper[i]<1.0e30)
     4247          upper[i] *= scale;
     4248      }
     4249    }
     4250    // scale partial matrix
     4251    {
     4252      CoinPackedMatrix * matrix = model2->matrix();
     4253      const int * row = matrix->getIndices();
     4254      const int * columnLength = matrix->getVectorLengths();
     4255      const CoinBigIndex * columnStart = matrix->getVectorStarts();
     4256      double * element = matrix->getMutableElements();
     4257      for (int i=0;i<numberNormal;i++) {
     4258        for (int j = columnStart[i];
     4259             j < columnStart[i] + columnLength[i]; j++) {
     4260          int iRow = row[j];
     4261          iRow = whichRows[iRow];
     4262          double scaleBy = scaleArray[iRow];
     4263          element[j] *= scaleBy;
     4264        }
     4265      }
     4266    }
     4267    // adjust rhs
     4268    double * rowLower = model2->rowLower();
     4269    double * rowUpper = model2->rowUpper();
     4270    for (int i=0;i<numberNonGub;i++) {
     4271      int iRow = whichRows[i];
     4272      rowLower[i] = lower[iRow];
     4273      rowUpper[i] = upper[iRow];
     4274    }
     4275    int numberGubColumns = numberColumns - put - numberFixed;
     4276    CoinBigIndex numberElements=0;
     4277    int * temp1 = new int [numberRows+1];
     4278    // get counts
     4279    memset(temp1,0,numberRows*sizeof(int));
     4280    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     4281      int iGub = columnIsGub[iColumn];
     4282      if (iGub>=0) {
     4283        numberElements += columnLength[iColumn]-1;
     4284        temp1[iGub]++;
     4285      }
     4286    }
     4287    /* Optional but means can eventually simplify coding
     4288       could even add in fixed slacks to deal with
     4289       singularities - but should not be necessary */
     4290    int numberSlacks=0;
     4291    for (int i = 0; i < numberRows; i++) {
     4292      if (rowIsGub[i]>=0) {
     4293        if (lower[i]<upper[i]) {
     4294          numberSlacks++;
     4295          temp1[i]++;
     4296        }
     4297      }
     4298    }
     4299    int * gubStart = new int [numberGub+1];
     4300    numberGub=0;
     4301    gubStart[0]=0;
     4302    for (int i = 0; i < numberRows; i++) {
     4303      if (rowIsGub[i]>=0) {
     4304        rowIsGub[i]=numberGub;
     4305        gubStart[numberGub+1]=gubStart[numberGub]+temp1[i];
     4306        temp1[numberGub]=0;
     4307        lower[numberGub]=lower[i];
     4308        upper[numberGub]=upper[i];
     4309        whichRows[numberNonGub+numberGub]=i;
     4310        numberGub++;
     4311      }
     4312    }
     4313    int numberGubColumnsPlus = numberGubColumns + numberSlacks;
     4314    double * lowerColumn2 = new double [numberGubColumnsPlus];
     4315    CoinFillN(lowerColumn2, numberGubColumnsPlus, 0.0);
     4316    double * upperColumn2 = new double [numberGubColumnsPlus];
     4317    CoinFillN(upperColumn2, numberGubColumnsPlus, COIN_DBL_MAX);
     4318    int * start2 = new int[numberGubColumnsPlus+1];
     4319    int * row2 = new int[numberElements];
     4320    double * element2 = new double[numberElements];
     4321    double * cost2 = new double [numberGubColumnsPlus];
     4322    CoinFillN(cost2, numberGubColumnsPlus, 0.0);
     4323    const double * cost = this->objective();
     4324    put = numberNormal;
     4325    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     4326      int iGub = columnIsGub[iColumn];
     4327      if (iGub>=0) {
     4328        // TEMP
     4329        //this->setColUpper(iColumn,COIN_DBL_MAX);
     4330        iGub = rowIsGub[iGub];
     4331        assert (iGub>=0);
     4332        int kPut = put+gubStart[iGub]+temp1[iGub];
     4333        temp1[iGub]++;
     4334        whichColumns[kPut]=iColumn;
     4335      }
     4336    }
     4337    for (int i = 0; i < numberRows; i++) {
     4338      if (rowIsGub[i]>=0) {
     4339        int iGub = rowIsGub[i];
     4340        if (lower[iGub]<upper[iGub]) {
     4341          int kPut = put+gubStart[iGub]+temp1[iGub];
     4342          temp1[iGub]++;
     4343          whichColumns[kPut]=iGub+numberColumns;
     4344        }
     4345      }
     4346    }
     4347    //this->primal(1); // TEMP
     4348    // redo rowIsGub to give lookup
     4349    for (int i=0;i<numberRows;i++)
     4350      rowIsGub[i]=-1;
     4351    for (int i=0;i<numberNonGub;i++)
     4352      rowIsGub[whichRows[i]]=i;
     4353    start2[0]=0;
     4354    numberElements = 0;
     4355    for (int i=0;i<numberGubColumnsPlus;i++) {
     4356      int iColumn = whichColumns[put++];
     4357      if (iColumn<numberColumns) {
     4358        cost2[i] = cost[iColumn];
     4359        lowerColumn2[i] = columnLower[iColumn];
     4360        upperColumn2[i] = columnUpper[iColumn];
     4361        upperColumn2[i] = COIN_DBL_MAX;
     4362        for (int j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     4363          int iRow = row[j];
     4364          double scaleBy = scaleArray[iRow];
     4365          iRow = rowIsGub[iRow];
     4366          if (iRow >= 0) {
     4367            row2[numberElements] = iRow;
     4368            element2[numberElements++] = elementByColumn[j]*scaleBy;
     4369          }
     4370        }
     4371      } else {
     4372        // slack
     4373        int iGub = iColumn-numberColumns;
     4374        double slack = upper[iGub]-lower[iGub];
     4375        assert (upper[iGub]<1.0e20);
     4376        lower[iGub]=upper[iGub];
     4377        cost2[i] = 0;
     4378        lowerColumn2[i] = 0;
     4379        upperColumn2[i] = slack;
     4380        upperColumn2[i] = COIN_DBL_MAX;
     4381      }
     4382      start2[i+1] = numberElements;
     4383    }
     4384    // clean up bounds on variables
     4385    for (int iSet=0;iSet<numberGub;iSet++) {
     4386      double lowerValue=0.0;
     4387      for (int i=gubStart[iSet];i<gubStart[iSet+1];i++) {
     4388        lowerValue += lowerColumn2[i];
     4389      }
     4390      assert (lowerValue<upper[iSet]+1.0e-6);
     4391      double gap = CoinMax(0.0,upper[iSet]-lowerValue);
     4392      for (int i=gubStart[iSet];i<gubStart[iSet+1];i++) {
     4393        if (upperColumn2[i]<1.0e30) {
     4394          upperColumn2[i] = CoinMin(upperColumn2[i],
     4395                                    lowerColumn2[i]+gap);
     4396        }
     4397      }
     4398    }
     4399    sprintf(message,"** Before adding matrix there are %d rows and %d columns",
     4400           model2->numberRows(), model2->numberColumns());
     4401    handler_->message(CLP_GENERAL, messages_)
     4402      << message << CoinMessageEol;
     4403    delete [] scaleArray;
     4404    delete [] temp1;
     4405    model2->setFactorizationFrequency(factorizationFrequency);
     4406    ClpDynamicMatrix * newMatrix =
     4407      new ClpDynamicMatrix(model2, numberGub,
     4408                                  numberGubColumnsPlus, gubStart,
     4409                                  lower, upper,
     4410                                  start2, row2, element2, cost2,
     4411                                  lowerColumn2, upperColumn2);
     4412    delete [] gubStart;
     4413    delete [] lowerColumn2;
     4414    delete [] upperColumn2;
     4415    delete [] start2;
     4416    delete [] row2;
     4417    delete [] element2;
     4418    delete [] cost2;
     4419    delete [] lower;
     4420    delete [] upper;
     4421    model2->replaceMatrix(newMatrix,true);
     4422#ifdef EVERY_ITERATION
     4423    {
     4424      ClpDynamicMatrix * gubMatrix =
     4425        dynamic_cast< ClpDynamicMatrix*>(model2->clpMatrix());
     4426      assert(gubMatrix);
     4427      gubMatrix->writeMps("gub.mps");
     4428    }
     4429#endif
     4430    delete [] columnIsGub;
     4431    delete [] rowIsGub;
     4432    newMatrix->switchOffCheck();
     4433#ifdef EVERY_ITERATION
     4434    newMatrix->setRefreshFrequency(1/*000*/);
     4435#else
     4436    newMatrix->setRefreshFrequency(1000);
     4437#endif
     4438    sprintf(message,
     4439            "** While after adding matrix there are %d rows and %d columns",
     4440            model2->numberRows(), model2->numberColumns());
     4441    handler_->message(CLP_GENERAL, messages_)
     4442      << message << CoinMessageEol;
     4443    model2->setSpecialOptions(4);    // exactly to bound
     4444    // Scaling off (done by hand)
     4445    model2->scaling(0);
     4446    return model2;
     4447  } else {
     4448    delete [] columnIsGub;
     4449    delete [] rowIsGub;
     4450    return NULL;
     4451  }
     4452}
     4453// Sets basis from original
     4454void
     4455ClpSimplexOther::setGubBasis(const ClpSimplex &original,const int * whichRows,
     4456                             const int * whichColumns)
     4457{
     4458  ClpDynamicMatrix * gubMatrix =
     4459    dynamic_cast< ClpDynamicMatrix*>(clpMatrix());
     4460  assert(gubMatrix);
     4461  int numberGubColumns = gubMatrix->numberGubColumns();
     4462  int numberNormal = gubMatrix->firstDynamic();
     4463  //int lastOdd = gubMatrix->firstAvailable();
     4464  //int numberTotalColumns = numberNormal + numberGubColumns;
     4465  //assert (numberTotalColumns==numberColumns+numberSlacks);
     4466  int numberRows = original.numberRows();
     4467  int numberColumns = original.numberColumns();
     4468  int * columnIsGub = new int [numberColumns];
     4469  int numberNonGub = gubMatrix->numberStaticRows();
     4470  //assert (firstOdd==numberNormal);
     4471  double * solution = primalColumnSolution();
     4472  const double * originalSolution = original.primalColumnSolution();
     4473  int numberSets = gubMatrix->numberSets();
     4474  const int * startSet = gubMatrix->startSets();
     4475  const CoinBigIndex * startColumn = gubMatrix->startColumn();
     4476  const double * columnLower = gubMatrix->columnLower();
     4477  for (int i=0;i<numberSets;i++) {
     4478    for (int j=startSet[i];j<startSet[i+1];j++) {
     4479      int iColumn = whichColumns[j+numberNormal];
     4480      if (iColumn<numberColumns)
     4481        columnIsGub[iColumn] = whichRows[numberNonGub+i];
     4482    }
     4483  }
     4484  int * numberKey = new int [numberRows];
     4485  memset(numberKey,0,numberRows*sizeof(int));
     4486  for (int i=0;i<numberGubColumns;i++) {
     4487    int iOrig = whichColumns[i+numberNormal];
     4488    if (iOrig<numberColumns) {
     4489      if (original.getColumnStatus(iOrig)==ClpSimplex::basic) {
     4490        int iRow = columnIsGub[iOrig];
     4491        assert (iRow>=0);
     4492        numberKey[iRow]++;
     4493      }
     4494    } else {
     4495      // Set slack
     4496      int iSet = iOrig - numberColumns;
     4497      int iRow = whichRows[iSet+numberNonGub];
     4498      if (original.getRowStatus(iRow)==ClpSimplex::basic)
     4499        numberKey[iRow]++;
     4500    }
     4501  }
     4502  /* Before going into cleanMatrix we need
     4503     gub status set (inSmall just means basic and active)
     4504     row status set
     4505  */
     4506  for (int i = 0; i < numberSets; i++) {
     4507    gubMatrix->setStatus(i,ClpSimplex::isFixed);
     4508  }
     4509  for (int i = 0; i < numberGubColumns; i++) {
     4510    int iOrig = whichColumns[i+numberNormal];
     4511    if (iOrig<numberColumns) {
     4512      ClpSimplex::Status status = original.getColumnStatus(iOrig);
     4513      if (status==ClpSimplex::atUpperBound) {
     4514        gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::atUpperBound);
     4515      } else if (status==ClpSimplex::atLowerBound) {
     4516        gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::atLowerBound);
     4517      } else if (status==ClpSimplex::basic) {
     4518        int iRow = columnIsGub[iOrig];
     4519        assert (iRow>=0);
     4520        assert(numberKey[iRow]);
     4521        if (numberKey[iRow]==1)
     4522          gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::soloKey);
     4523        else
     4524          gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::inSmall);
     4525      }
     4526    } else {
     4527      // slack
     4528      int iSet = iOrig - numberColumns;
     4529      int iRow = whichRows[iSet+numberNonGub];
     4530      if (original.getRowStatus(iRow)==ClpSimplex::basic) {
     4531        assert(numberKey[iRow]);
     4532        if (numberKey[iRow]==1)
     4533          gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::soloKey);
     4534        else
     4535          gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::inSmall);
     4536      } else {
     4537        gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::atLowerBound);
     4538      }
     4539    }
     4540  }
     4541  // deal with sets without key
     4542  for (int i = 0; i < numberSets; i++) {
     4543    int iRow = whichRows[numberNonGub+i];
     4544    if (!numberKey[iRow]) {
     4545      if (original.getRowStatus(iRow)==ClpSimplex::basic) {
     4546        gubMatrix->setStatus(i,ClpSimplex::basic);
     4547      } else {
     4548        // If not at lb make key otherwise one with smallest number els
     4549        double largest=0.0;
     4550        int fewest=numberRows+1;
     4551        int chosen=-1;
     4552        for (int j=startSet[i];j<startSet[i+1];j++) {
     4553          int length=startColumn[j+1]-startColumn[j];
     4554          int iOrig = whichColumns[j+numberNormal];
     4555          double value;
     4556          if (iOrig<numberColumns) {
     4557            value=originalSolution[iOrig]-columnLower[j];
     4558          } else {
     4559            // slack - take value as 0.0 as will win on length
     4560            value=0.0;
     4561          }
     4562          if (value>largest+1.0e-8) {
     4563            largest=value;
     4564            fewest=length;
     4565            chosen=j;
     4566          } else if (fabs(value-largest)<=1.0e-8&&length<fewest) {
     4567            largest=value;
     4568            fewest=length;
     4569            chosen=j;
     4570          }
     4571        }
     4572        assert(chosen>=0);
     4573        // set as key
     4574        for (int j=startSet[i];j<startSet[i+1];j++) {
     4575          if (j!=chosen)
     4576            gubMatrix->setDynamicStatus(j,ClpDynamicMatrix::atLowerBound);
     4577          else
     4578            gubMatrix->setDynamicStatus(j,ClpDynamicMatrix::soloKey);
     4579        }
     4580      }
     4581    }
     4582  }
     4583  for (int i = 0; i < numberNormal; i++) {
     4584    int iOrig = whichColumns[i];
     4585    setColumnStatus(i,original.getColumnStatus(iOrig));
     4586    solution[i]=originalSolution[iOrig];
     4587  }
     4588  for (int i = 0; i < numberNonGub; i++) {
     4589    int iOrig = whichRows[i];
     4590    setRowStatus(i,original.getRowStatus(iOrig));
     4591  }
     4592  // Fill in current matrix
     4593  gubMatrix->initialProblem();
     4594  delete [] numberKey;
     4595  delete [] columnIsGub;
     4596}
     4597// Restores basis to original
     4598void
     4599ClpSimplexOther::getGubBasis(ClpSimplex &original,const int * whichRows,
     4600                             const int * whichColumns) const
     4601{
     4602  ClpDynamicMatrix * gubMatrix =
     4603    dynamic_cast< ClpDynamicMatrix*>(clpMatrix());
     4604  assert(gubMatrix);
     4605  int numberGubColumns = gubMatrix->numberGubColumns();
     4606  int numberNormal = gubMatrix->firstDynamic();
     4607  //int lastOdd = gubMatrix->firstAvailable();
     4608  //int numberRows = original.numberRows();
     4609  int numberColumns = original.numberColumns();
     4610  int numberNonGub = gubMatrix->numberStaticRows();
     4611  //assert (firstOdd==numberNormal);
     4612  double * solution = primalColumnSolution();
     4613  double * originalSolution = original.primalColumnSolution();
     4614  int numberSets = gubMatrix->numberSets();
     4615  const double * cost = original.objective();
     4616  int lastOdd = gubMatrix->firstAvailable();
     4617  //assert (numberTotalColumns==numberColumns+numberSlacks);
     4618  int numberRows = original.numberRows();
     4619  //int numberStaticRows = gubMatrix->numberStaticRows();
     4620  const int * startSet = gubMatrix->startSets();
     4621  unsigned char * status = original.statusArray();
     4622  unsigned char * rowStatus = status+numberColumns;
     4623  //assert (firstOdd==numberNormal);
     4624  for (int i=0;i<numberSets;i++) {
     4625    int iRow = whichRows[i+numberNonGub];
     4626    original.setRowStatus(iRow,ClpSimplex::atLowerBound);
     4627  }
     4628  const int * id = gubMatrix->id();
     4629  const double * columnLower = gubMatrix->columnLower();
     4630  const double * columnUpper = gubMatrix->columnUpper();
     4631  for (int i = 0; i < numberGubColumns; i++) {
     4632    int iOrig = whichColumns[i+numberNormal];
     4633    if (iOrig<numberColumns) {
     4634      if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::atUpperBound) {
     4635        originalSolution[iOrig] = columnUpper[i];
     4636        status[iOrig] = 2;
     4637      } else if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::atLowerBound && columnLower) {
     4638        originalSolution[iOrig] = columnLower[i];
     4639        status[iOrig] = 3;
     4640      } else if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::soloKey) {
     4641        int iSet = gubMatrix->whichSet(i);
     4642        originalSolution[iOrig] = gubMatrix->keyValue(iSet);
     4643        status[iOrig] = 1;
     4644      } else {
     4645        originalSolution[iOrig] = 0.0;
     4646        status[iOrig] = 4;
     4647      }
     4648    } else {
     4649      // slack
     4650      int iSet = iOrig - numberColumns;
     4651      int iRow = whichRows[iSet+numberNonGub];
     4652      if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::atUpperBound) {
     4653        original.setRowStatus(iRow,ClpSimplex::atLowerBound);
     4654      } else if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::atLowerBound) {
     4655        original.setRowStatus(iRow,ClpSimplex::atUpperBound);
     4656      } else if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::soloKey) {
     4657        original.setRowStatus(iRow,ClpSimplex::basic);
     4658      }
     4659    }
     4660  }
     4661  for (int i = 0; i < numberNormal; i++) {
     4662    int iOrig = whichColumns[i];
     4663    ClpSimplex::Status thisStatus = getStatus(i);
     4664    if (thisStatus == ClpSimplex::basic)
     4665      status[iOrig] = 1;
     4666    else if (thisStatus == ClpSimplex::atLowerBound)
     4667      status[iOrig] = 3;
     4668    else if (thisStatus == ClpSimplex::atUpperBound)
     4669      status[iOrig] = 2;
     4670    else if (thisStatus == ClpSimplex::isFixed)
     4671      status[iOrig] = 5;
     4672    else
     4673      abort();
     4674    originalSolution[iOrig] = solution[i];
     4675  }
     4676  for (int i = numberNormal; i < lastOdd; i++) {
     4677    int iOrig = whichColumns[id[i-numberNormal] + numberNormal];
     4678    if (iOrig<numberColumns) {
     4679      ClpSimplex::Status thisStatus = getStatus(i);
     4680      if (thisStatus == ClpSimplex::basic)
     4681        status[iOrig] = 1;
     4682      else if (thisStatus == ClpSimplex::atLowerBound)
     4683        status[iOrig] = 3;
     4684      else if (thisStatus == ClpSimplex::atUpperBound)
     4685        status[iOrig] = 2;
     4686      else if (thisStatus == ClpSimplex::isFixed)
     4687        status[iOrig] = 5;
     4688      else
     4689        abort();
     4690      originalSolution[iOrig] = solution[i];
     4691    } else {
     4692      // slack (basic probably)
     4693      int iSet = iOrig - numberColumns;
     4694      int iRow = whichRows[iSet+numberNonGub];
     4695      ClpSimplex::Status thisStatus = getStatus(i);
     4696      if (thisStatus == ClpSimplex::atLowerBound)
     4697        thisStatus = ClpSimplex::atUpperBound;
     4698      else if (thisStatus == ClpSimplex::atUpperBound)
     4699        thisStatus = ClpSimplex::atLowerBound;
     4700      original.setRowStatus(iRow,thisStatus);
     4701    }
     4702  }
     4703  for (int i = 0; i < numberNonGub; i++) {
     4704    int iOrig = whichRows[i];
     4705    ClpSimplex::Status thisStatus = getRowStatus(i);
     4706    if (thisStatus == ClpSimplex::basic)
     4707      rowStatus[iOrig] = 1;
     4708    else if (thisStatus == ClpSimplex::atLowerBound)
     4709      rowStatus[iOrig] = 3;
     4710    else if (thisStatus == ClpSimplex::atUpperBound)
     4711      rowStatus[iOrig] = 2;
     4712    else if (thisStatus == ClpSimplex::isFixed)
     4713      rowStatus[iOrig] = 5;
     4714    else
     4715      abort();
     4716  }
     4717  int * numberKey = new int [numberRows];
     4718  memset(numberKey,0,numberRows*sizeof(int));
     4719  for (int i=0;i<numberSets;i++) {
     4720    int iRow = whichRows[i+numberNonGub];
     4721    for (int j=startSet[i];j<startSet[i+1];j++) {
     4722      int iOrig = whichColumns[j+numberNormal];
     4723      if (iOrig<numberColumns) {
     4724        if (original.getColumnStatus(iOrig)==ClpSimplex::basic) {
     4725          numberKey[iRow]++;
     4726        }
     4727      } else {
     4728        // slack
     4729        if (original.getRowStatus(iRow)==ClpSimplex::basic)
     4730          numberKey[iRow]++;
     4731      }
     4732    }
     4733  }
     4734  for (int i=0;i<numberSets;i++) {
     4735    int iRow = whichRows[i+numberNonGub];
     4736    if (!numberKey[iRow]) {
     4737      original.setRowStatus(iRow,ClpSimplex::basic);
     4738    }
     4739  }
     4740  delete [] numberKey;
     4741  double objValue = 0.0;
     4742  for (int i = 0; i < numberColumns; i++)
     4743    objValue += cost[i] * originalSolution[i];
     4744  //printf("objective value is %g\n", objValue);
     4745}
  • trunk/Clp/src/ClpSimplexOther.hpp

    r1665 r1676  
    190190                      const int * whichRows, const int * whichColumns,
    191191                      int nBound);
     192     /** Returns gub version of model or NULL
     193         whichRows has to be numberRows
     194         whichColumns has to be numberRows+numberColumns */
     195     ClpSimplex * gubVersion(int * whichRows, int * whichColumns,
     196                             int neededGub,
     197                             int factorizationFrequency=50);
     198     /// Sets basis from original
     199     void setGubBasis(const ClpSimplex &original,const int * whichRows,
     200                      const int * whichColumns);
     201     /// Restores basis to original
     202     void getGubBasis(ClpSimplex &original,const int * whichRows,
     203                      const int * whichColumns) const;
    192204     /// Quick try at cleaning up duals if postsolve gets wrong
    193205     void cleanupAfterPostsolve();
  • trunk/Clp/src/ClpSimplexPrimal.cpp

    r1671 r1676  
    10921092          objectiveWork_ = cost_;
    10931093          rowObjectiveWork_ = cost_ + numberColumns_;
    1094           if (n)
     1094          if (n||matrix_->type()>=15)
    10951095               gutsOfSolution(NULL, NULL, ifValuesPass != 0);
    10961096     }
     
    11181118                    && firstFree_ < 0) {
    11191119               if (handler_->logLevel() == 63)
    1120                     printf("lastobj %g this %g force %d ", lastObj, thisObj, forceFactorization_);
     1120                    printf("lastobj %g this %g force %d\n", lastObj, thisObj, forceFactorization_);
    11211121               int maxFactor = factorization_->maximumPivots();
    11221122               if (maxFactor > 10) {
     
    11301130                     && firstFree_ < 0) {
    11311131               if (handler_->logLevel() == 63)
    1132                     printf("lastobj3 %g this3 %g `force %d ", lastObj3, thisObj, forceFactorization_);
     1132                    printf("lastobj3 %g this3 %g force %d\n", lastObj3, thisObj, forceFactorization_);
    11331133               int maxFactor = factorization_->maximumPivots();
    11341134               if (maxFactor > 10) {
  • trunk/Clp/src/ClpSolve.cpp

    r1665 r1676  
    30253025     if (matched == (1 << (CLP_PROGRESS - 1)))
    30263026          numberMatched = 0;
    3027      if (numberMatched) {
     3027     if (numberMatched && model_->clpMatrix()->type() < 15) {
    30283028          model_->messageHandler()->message(CLP_POSSIBLELOOP, model_->messages())
    30293029                    << numberMatched
Note: See TracChangeset for help on using the changeset viewer.