Changeset 1330 for trunk/Clp


Ignore:
Timestamp:
Feb 6, 2009 8:55:11 AM (11 years ago)
Author:
forrest
Message:

improve dualize

File:
1 edited

Legend:

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

    r1321 r1330  
    865865  }
    866866  int kRow=0;
     867  int kExtraRow=numberRows;
    867868  for (iRow=0;iRow<numberRows;iRow++) {
    868869    if (rowLower[iRow]<-1.0e20) {
     
    893894        which[kRow]=iRow;
    894895        kRow++;
    895         newObjective[kRow]=-rowLower[iRow];
    896         fromRowsLower[kRow]=0.0;
    897         fromRowsUpper[kRow]=COIN_DBL_MAX;
    898         which[kRow]=iRow;
    899         kRow++;
     896        newObjective[kExtraRow]=-rowLower[iRow];
     897        fromRowsLower[kExtraRow]=0.0;
     898        fromRowsUpper[kExtraRow]=COIN_DBL_MAX;
     899        which[kExtraRow]=iRow;
     900        kExtraRow++;
    900901      }
    901902    }
     
    905906    newCopy.setExtraGap(0.0);
    906907    newCopy.setExtraMajor(0.0);
    907     newCopy.submatrixOfWithDuplicates(rowCopy,kRow,which);
     908    newCopy.submatrixOfWithDuplicates(rowCopy,kExtraRow,which);
    908909    rowCopy=newCopy;
    909910  }
     
    959960  const double * dualActs = dualProblem->primalRowSolution();
    960961#if 0
    961   const double * primalDual = this->dualRowSolution();
    962   const double * primalDj = this->dualColumnSolution();
    963   const double * primalSol = this->primalColumnSolution();
    964   const double * primalActs = this->primalRowSolution();
     962  ClpSimplex thisCopy=*this;
     963  thisCopy.dual(); // for testing
     964  const double * primalDual = thisCopy.dualRowSolution();
     965  const double * primalDj = thisCopy.dualColumnSolution();
     966  const double * primalSol = thisCopy.primalColumnSolution();
     967  const double * primalActs = thisCopy.primalRowSolution();
    965968  char ss[]={'F','B','U','L','S','F'};
    966   dual(); // for testing
    967969  printf ("Dual problem row info %d rows\n",dualProblem->numberRows());
    968970  for (iRow=0;iRow<dualProblem->numberRows();iRow++)
     
    975977           iColumn,ss[dualProblem->getColumnStatus(iColumn)],
    976978           dualSol[iColumn],dualDj[iColumn]);
    977   printf ("Primal problem row info %d rows\n",this->numberRows());
    978   for (iRow=0;iRow<this->numberRows();iRow++)
     979  printf ("Primal problem row info %d rows\n",thisCopy.numberRows());
     980  for (iRow=0;iRow<thisCopy.numberRows();iRow++)
    979981    printf("%d at %c primal %g dual %g\n",
    980            iRow,ss[this->getRowStatus(iRow)],
     982           iRow,ss[thisCopy.getRowStatus(iRow)],
    981983           primalActs[iRow],primalDual[iRow]);
    982   printf ("Primal problem column info %d columns\n",this->numberColumns());
    983   for (iColumn=0;iColumn<this->numberColumns();iColumn++)
     984  printf ("Primal problem column info %d columns\n",thisCopy.numberColumns());
     985  for (iColumn=0;iColumn<thisCopy.numberColumns();iColumn++)
    984986    printf("%d at %c primal %g dual %g\n",
    985            iColumn,ss[this->getColumnStatus(iColumn)],
     987           iColumn,ss[thisCopy.getColumnStatus(iColumn)],
    986988           primalSol[iColumn],primalDj[iColumn]);
    987989#endif
    988990  // position at bound information
    989   int jColumn=numberRows_+numberExtraRows;
     991  int jColumn=numberRows_;
    990992  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    991993    double objValue =optimizationDirection_*objective[iColumn];
     
    10341036        } else if (columnUpper_[iColumn]<1.0e20) {
    10351037          columnActivity_[iColumn]=-dualDual[iColumn] + columnUpper_[iColumn];
     1038        } else {
     1039          columnActivity_[iColumn]=-dualDual[iColumn];
    10361040        }
    10371041        reducedCost_[iColumn]=0.0;
     
    10431047          setColumnStatus(iColumn,basic);
    10441048          numberBasic++;
     1049          //printf("Col %d otherV %g dualDual %g\n",iColumn,
     1050          // otherValue,dualDual[iColumn]);
    10451051          columnActivity_[iColumn]=-dualDual[iColumn];
     1052          columnActivity_[iColumn]=otherValue;
    10461053          reducedCost_[iColumn]=0.0;
    10471054        } else {
     
    10621069  // now rows
    10631070  int kRow=0;
     1071  int kExtraRow=jColumn;
    10641072  int numberRanges=0;
    10651073  for (iRow=0;iRow<numberRows_;iRow++) {
     
    11001108        // range
    11011109        numberRanges++;
    1102         Status statusL = dualProblem->getColumnStatus(kRow+1);
     1110        Status statusL = dualProblem->getColumnStatus(kExtraRow);
     1111        kExtraRow++;
    11031112        if (status==basic) {
    11041113          assert (statusL!=basic);
     
    11161125        }
    11171126        kRow++;
    1118         kRow++;
    1119       }
    1120     }
    1121   }
    1122   if (numberRanges) {
    1123     printf("%d ranges - coding needed\n",numberRanges);
     1127      }
     1128    }
     1129  }
     1130  if (numberBasic!=numberRows_) {
     1131    printf("Bad basis - ranges - coding needed\n");
     1132    assert (numberRanges);
    11241133    returnCode=1;
    1125   }
    1126   if (numberBasic!=numberRows_) {
    1127     printf("Bad basis - ranges?\n");
    1128     assert (numberRanges);
    11291134  }
    11301135  if (optimizationDirection_<0.0) {
     
    11321137      dual_[iRow]=-dual_[iRow];
    11331138    }
    1134     for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    1135       reducedCost_[iColumn]=-reducedCost_[iColumn];
    1136     }
    11371139  }
    11381140  // redo row activities
    11391141  memset(rowActivity_,0,numberRows_*sizeof(double));
    1140   matrix_->times(-1.0,columnActivity_,rowActivity_);
     1142  matrix_->times(1.0,columnActivity_,rowActivity_);
     1143  // redo reduced costs
     1144  memcpy(reducedCost_,this->objective(),numberColumns_*sizeof(double));
     1145  matrix_->transposeTimes(-1.0,dual_,reducedCost_);
    11411146  checkSolutionInternal();
    1142   return 1; //temp
     1147  // Below will go to ..DEBUG later
     1148  if (returnCode)
     1149    return 1; // Skip checks as will fail at present
     1150#ifndef NDEBUG
     1151  // Check if correct
     1152  double * columnActivity = CoinCopyOfArray(columnActivity_,numberColumns_);
     1153  double * rowActivity = CoinCopyOfArray(rowActivity_,numberRows_);
     1154  double * reducedCost = CoinCopyOfArray(reducedCost_,numberColumns_);
     1155  double * dual = CoinCopyOfArray(dual_,numberRows_);
     1156  this->dual();
     1157  CoinRelFltEq eq(1.0e-5);
     1158  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     1159    assert(eq(columnActivity[iColumn],columnActivity_[iColumn]));
     1160  }
     1161  for (iRow=0;iRow<numberRows_;iRow++) {
     1162    assert(eq(dual[iRow],dual_[iRow]));
     1163  }
     1164  for (iRow=0;iRow<numberRows_;iRow++) {
     1165    assert(eq(rowActivity[iRow],rowActivity_[iRow]));
     1166  }
     1167  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     1168    assert(eq(reducedCost[iColumn],reducedCost_[iColumn]));
     1169  }
     1170  delete [] columnActivity;
     1171  delete [] rowActivity;
     1172  delete [] reducedCost;
     1173  delete [] dual;
     1174#endif
    11431175  return returnCode;
    11441176}
Note: See TracChangeset for help on using the changeset viewer.