Changeset 1147


Ignore:
Timestamp:
Jan 2, 2008 5:14:29 AM (12 years ago)
Author:
forrest
Message:

changes to try and make faster

Location:
trunk/Clp/src
Files:
18 edited

Legend:

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

    r1140 r1147  
    16371637    parameters[numberParameters++]=
    16381638      CbcOrClpParam("force!Solution","Whether to use given solution as crash for BAB",
    1639                     "off",USESOLUTION);
    1640     parameters[numberParameters-1].append("on");
    1641   parameters[numberParameters-1].setLonghelp
    1642     (
    1643      "If on then tries to branch to solution given by AMPL or priorities file."
     1639                    -1,20000000,USESOLUTION);
     1640    parameters[numberParameters-1].setIntValue(-1);
     1641    parameters[numberParameters-1].setLonghelp
     1642    (
     1643     "-1 off.  If 0 then tries to branch to solution given by AMPL or priorities file. \
     1644If >0 then also does that many nodes on fixed problem."
    16441645     );
    16451646#endif
  • trunk/Clp/src/CbcOrClpParam.hpp

    r1137 r1147  
    6161    MAXFACTOR,PERTVALUE,MAXITERATION,PRESOLVEPASS,IDIOT,SPRINT,
    6262    OUTPUTFORMAT,SLPVALUE,PRESOLVEOPTIONS,PRINTOPTIONS,SPECIALOPTIONS,
    63     SUBSTITUTION,DUALIZE,VERBOSE,CPP,PROCESSTUNE,
     63    SUBSTITUTION,DUALIZE,VERBOSE,CPP,PROCESSTUNE,USESOLUTION,
    6464
    6565    STRONGBRANCHING=151,CUTDEPTH, MAXNODES,NUMBERBEFORE,NUMBERANALYZE,
     
    7878    GOMORYCUTS,PROBINGCUTS,KNAPSACKCUTS,REDSPLITCUTS,
    7979    ROUNDING,SOLVER,CLIQUECUTS,COSTSTRATEGY,FLOWCUTS,MIXEDCUTS,
    80     TWOMIRCUTS,PREPROCESS,FPUMP,GREEDY,COMBINE,LOCALTREE,USESOLUTION,SOS,
     80    TWOMIRCUTS,PREPROCESS,FPUMP,GREEDY,COMBINE,LOCALTREE,SOS,
    8181    LANDPCUTS,RINS,RESIDCUTS,RENS,
    8282   
  • trunk/Clp/src/ClpDualRowSteepest.cpp

    r1141 r1147  
    388388    const int *permute = model_->factorization()->permute();
    389389    //double * region = alternateWeights_->denseVector();
    390     for ( i = 0; i < numberNonZero; i ++ ) {
    391       int iRow = which[i];
    392       double value = work[i];
    393       norm += value*value;
    394       iRow = permute[iRow];
    395       work2[iRow] = value;
    396       which2[i] = iRow;
     390    if (permute) {
     391      for ( i = 0; i < numberNonZero; i ++ ) {
     392        int iRow = which[i];
     393        double value = work[i];
     394        norm += value*value;
     395        iRow = permute[iRow];
     396        work2[iRow] = value;
     397        which2[i] = iRow;
     398      }
     399    } else {
     400      for ( i = 0; i < numberNonZero; i ++ ) {
     401        int iRow = which[i];
     402        double value = work[i];
     403        norm += value*value;
     404        //iRow = permute[iRow];
     405        work2[iRow] = value;
     406        which2[i] = iRow;
     407      }
    397408    }
    398409    spare->setNumElements ( numberNonZero );
  • trunk/Clp/src/ClpFactorization.cpp

    r1111 r1147  
    2424// Constructors / Destructor / Assignment
    2525//#############################################################################
     26#ifndef CLP_MULTIPLE_FACTORIZATIONS   
    2627
    2728//-------------------------------------------------------------------
     
    997998  }
    998999}
     1000#else
     1001// This one allows multiple factorizations
     1002
     1003//-------------------------------------------------------------------
     1004// Default Constructor
     1005//-------------------------------------------------------------------
     1006ClpFactorization::ClpFactorization ()
     1007{
     1008#ifndef SLIM_CLP
     1009  networkBasis_ = NULL;
     1010#endif
     1011  coinFactorizationA_ = new CoinFactorization() ;
     1012}
     1013
     1014//-------------------------------------------------------------------
     1015// Copy constructor
     1016//-------------------------------------------------------------------
     1017ClpFactorization::ClpFactorization (const ClpFactorization & rhs)
     1018{
     1019#ifndef SLIM_CLP
     1020  if (rhs.networkBasis_)
     1021    networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
     1022  else
     1023    networkBasis_=NULL;
     1024#endif
     1025  if (rhs.coinFactorizationA_)
     1026    coinFactorizationA_ = new CoinFactorization(*(rhs.coinFactorizationA_));
     1027  else
     1028    coinFactorizationA_=NULL;
     1029}
     1030
     1031ClpFactorization::ClpFactorization (const CoinFactorization & rhs)
     1032{
     1033#ifndef SLIM_CLP
     1034  networkBasis_=NULL;
     1035#endif
     1036  coinFactorizationA_ = new CoinFactorization(rhs);
     1037}
     1038
     1039//-------------------------------------------------------------------
     1040// Destructor
     1041//-------------------------------------------------------------------
     1042ClpFactorization::~ClpFactorization ()
     1043{
     1044#ifndef SLIM_CLP
     1045  delete networkBasis_;
     1046#endif
     1047  delete coinFactorizationA_;
     1048}
     1049
     1050//----------------------------------------------------------------
     1051// Assignment operator
     1052//-------------------------------------------------------------------
     1053ClpFactorization &
     1054ClpFactorization::operator=(const ClpFactorization& rhs)
     1055{
     1056  if (this != &rhs) {
     1057#ifndef SLIM_CLP
     1058    delete networkBasis_;
     1059    if (rhs.networkBasis_)
     1060      networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
     1061    else
     1062      networkBasis_=NULL;
     1063#endif
     1064    delete coinFactorizationA_;
     1065    if (rhs.coinFactorizationA_)
     1066      coinFactorizationA_ = new CoinFactorization(*(rhs.coinFactorizationA_));
     1067    else
     1068      coinFactorizationA_=NULL;
     1069  }
     1070  return *this;
     1071}
     1072#if 0
     1073static unsigned int saveList[10000];
     1074int numberSave=-1;
     1075inline bool isDense(int i) {
     1076  return ((saveList[i>>5]>>(i&31))&1)!=0;
     1077}
     1078inline void setDense(int i) {
     1079  unsigned int & value = saveList[i>>5];
     1080  int bit = i&31;
     1081  value |= (1<<bit);
     1082}
     1083#endif
     1084int
     1085ClpFactorization::factorize ( ClpSimplex * model,
     1086                              int solveType, bool valuesPass)
     1087{
     1088  ClpMatrixBase * matrix = model->clpMatrix();
     1089  int numberRows = model->numberRows();
     1090  int numberColumns = model->numberColumns();
     1091  if (!numberRows)
     1092    return 0;
     1093  // If too many compressions increase area
     1094  if (coinFactorizationA_->pivots()>1&&coinFactorizationA_->numberCompressions()*10 > coinFactorizationA_->pivots()+10) {
     1095    coinFactorizationA_->areaFactor( coinFactorizationA_->areaFactor()* 1.1);
     1096  }
     1097  //int numberPivots=coinFactorizationA_->pivots();
     1098#if 0
     1099  if (model->algorithm()>0)
     1100    numberSave=-1;
     1101#endif
     1102#ifndef SLIM_CLP
     1103  if (!networkBasis_||doCheck) {
     1104#endif
     1105    coinFactorizationA_->setStatus(-99);
     1106    int * pivotVariable = model->pivotVariable();
     1107    //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
     1108    while (coinFactorizationA_->status()<-98) {
     1109     
     1110      int i;
     1111      int numberBasic=0;
     1112      int numberRowBasic;
     1113      // Move pivot variables across if they look good
     1114      int * pivotTemp = model->rowArray(0)->getIndices();
     1115      assert (!model->rowArray(0)->getNumElements());
     1116      if (!matrix->rhsOffset(model)) {
     1117#if 0
     1118        if (numberSave>0) {
     1119          int nStill=0;
     1120          int nAtBound=0;
     1121          int nZeroDual=0;
     1122          CoinIndexedVector * array = model->rowArray(3);
     1123          CoinIndexedVector * objArray = model->columnArray(1);
     1124          array->clear();
     1125          objArray->clear();
     1126          double * cost = model->costRegion();
     1127          double tolerance = model->primalTolerance();
     1128          double offset=0.0;
     1129          for (i=0;i<numberRows;i++) {
     1130            int iPivot = pivotVariable[i];
     1131            if (iPivot<numberColumns&&isDense(iPivot)) {
     1132              if (model->getColumnStatus(iPivot)==ClpSimplex::basic) {
     1133                nStill++;
     1134                double value=model->solutionRegion()[iPivot];
     1135                double dual = model->dualRowSolution()[i];
     1136                double lower=model->lowerRegion()[iPivot];
     1137                double upper=model->upperRegion()[iPivot];
     1138                ClpSimplex::Status status;
     1139                if (fabs(value-lower)<tolerance) {
     1140                  status=ClpSimplex::atLowerBound;
     1141                  nAtBound++;
     1142                } else if (fabs(value-upper)<tolerance) {
     1143                  nAtBound++;
     1144                  status=ClpSimplex::atUpperBound;
     1145                } else if (value>lower&&value<upper) {
     1146                  status=ClpSimplex::superBasic;
     1147                } else {
     1148                  status=ClpSimplex::basic;
     1149                }
     1150                if (status!=ClpSimplex::basic) {
     1151                  if (model->getRowStatus(i)!=ClpSimplex::basic) {
     1152                    model->setColumnStatus(iPivot,ClpSimplex::atLowerBound);
     1153                    model->setRowStatus(i,ClpSimplex::basic);
     1154                    pivotVariable[i]=i+numberColumns;
     1155                    model->dualRowSolution()[i]=0.0;
     1156                    model->djRegion(0)[i]=0.0;
     1157                    array->add(i,dual);
     1158                    offset += dual*model->solutionRegion(0)[i];
     1159                  }
     1160                }
     1161                if (fabs(dual)<1.0e-5)
     1162                  nZeroDual++;
     1163              }
     1164            }
     1165          }
     1166          printf("out of %d dense, %d still in basis, %d at bound, %d with zero dual - offset %g\n",
     1167                 numberSave,nStill,nAtBound,nZeroDual,offset);
     1168          if (array->getNumElements()) {
     1169            // modify costs
     1170            model->clpMatrix()->transposeTimes(model,1.0,array,model->columnArray(0),
     1171                                               objArray);
     1172            array->clear();
     1173            int n=objArray->getNumElements();
     1174            int * indices = objArray->getIndices();
     1175            double * elements = objArray->denseVector();
     1176            for (i=0;i<n;i++) {
     1177              int iColumn = indices[i];
     1178              cost[iColumn] -= elements[iColumn];
     1179              elements[iColumn]=0.0;
     1180            }
     1181            objArray->setNumElements(0);
     1182          }
     1183        }
     1184#endif
     1185        // Seems to prefer things in order so quickest
     1186        // way is to go though like this
     1187        for (i=0;i<numberRows;i++) {
     1188          if (model->getRowStatus(i) == ClpSimplex::basic)
     1189            pivotTemp[numberBasic++]=i;
     1190        }
     1191        numberRowBasic=numberBasic;
     1192        /* Put column basic variables into pivotVariable
     1193           This is done by ClpMatrixBase to allow override for gub
     1194        */
     1195        matrix->generalExpanded(model,0,numberBasic);
     1196      } else {
     1197        // Long matrix - do a different way
     1198        bool fullSearch=false;
     1199        for (i=0;i<numberRows;i++) {
     1200          int iPivot = pivotVariable[i];
     1201          if (iPivot>=numberColumns) {
     1202            pivotTemp[numberBasic++]=iPivot-numberColumns;
     1203          }
     1204        }
     1205        numberRowBasic=numberBasic;
     1206        for (i=0;i<numberRows;i++) {
     1207          int iPivot = pivotVariable[i];
     1208          if (iPivot<numberColumns) {
     1209            if (iPivot>=0) {
     1210              pivotTemp[numberBasic++]=iPivot;
     1211            } else {
     1212              // not full basis
     1213              fullSearch=true;
     1214              break;
     1215            }
     1216          }
     1217        }
     1218        if (fullSearch) {
     1219          // do slow way
     1220          numberBasic=0;
     1221          for (i=0;i<numberRows;i++) {
     1222            if (model->getRowStatus(i) == ClpSimplex::basic)
     1223              pivotTemp[numberBasic++]=i;
     1224          }
     1225          numberRowBasic=numberBasic;
     1226          /* Put column basic variables into pivotVariable
     1227             This is done by ClpMatrixBase to allow override for gub
     1228          */
     1229          matrix->generalExpanded(model,0,numberBasic);
     1230        }
     1231      }
     1232      if (numberBasic>model->maximumBasic()) {
     1233#if 0 // ndef NDEBUG
     1234        printf("%d basic - should only be %d\n",
     1235               numberBasic,numberRows);
     1236#endif
     1237        // Take out some
     1238        numberBasic=numberRowBasic;
     1239        for (int i=0;i<numberColumns;i++) {
     1240          if (model->getColumnStatus(i) == ClpSimplex::basic) {
     1241            if (numberBasic<numberRows)
     1242              numberBasic++;
     1243            else
     1244              model->setColumnStatus(i,ClpSimplex::superBasic);
     1245          }
     1246        }
     1247        numberBasic=numberRowBasic;
     1248        matrix->generalExpanded(model,0,numberBasic);
     1249      }
     1250#ifndef SLIM_CLP
     1251      // see if matrix a network
     1252#ifndef NO_RTTI
     1253      ClpNetworkMatrix* networkMatrix =
     1254        dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
     1255#else
     1256      ClpNetworkMatrix* networkMatrix = NULL;
     1257      if (model->clpMatrix()->type()==11)
     1258        networkMatrix =
     1259        static_cast< ClpNetworkMatrix*>(model->clpMatrix());
     1260#endif
     1261      // If network - still allow ordinary factorization first time for laziness
     1262      if (networkMatrix)
     1263        coinFactorizationA_->setBiasLU(0); // All to U if network
     1264      //int saveMaximumPivots = maximumPivots();
     1265      delete networkBasis_;
     1266      networkBasis_ = NULL;
     1267      if (networkMatrix&&!doCheck)
     1268        maximumPivots(1);
     1269#endif
     1270      //printf("L, U, R %d %d %d\n",numberElementsL(),numberElementsU(),numberElementsR());
     1271      while (coinFactorizationA_->status()==-99) {
     1272        // maybe for speed will be better to leave as many regions as possible
     1273        coinFactorizationA_->gutsOfDestructor();
     1274        coinFactorizationA_->gutsOfInitialize(2);
     1275        CoinBigIndex numberElements=numberRowBasic;
     1276
     1277        // compute how much in basis
     1278
     1279        int i;
     1280        // can change for gub
     1281        int numberColumnBasic = numberBasic-numberRowBasic;
     1282
     1283        numberElements +=matrix->countBasis(model,
     1284                                           pivotTemp+numberRowBasic,
     1285                                           numberRowBasic,
     1286                                            numberColumnBasic);
     1287        // and recompute as network side say different
     1288        if (model->numberIterations())
     1289          numberRowBasic = numberBasic - numberColumnBasic;
     1290        numberElements = 3 * numberBasic + 3 * numberElements + 10000;
     1291        coinFactorizationA_->getAreas ( numberRows,
     1292                   numberRowBasic+numberColumnBasic, numberElements,
     1293                   2 * numberElements );
     1294        //fill
     1295        // Fill in counts so we can skip part of preProcess
     1296        int * numberInRow = coinFactorizationA_->numberInRow();
     1297        int * numberInColumn = coinFactorizationA_->numberInColumn();
     1298        CoinZeroN ( numberInRow, coinFactorizationA_->numberRows() + 1 );
     1299        CoinZeroN ( numberInColumn, coinFactorizationA_->maximumColumnsExtra() + 1 );
     1300        double * elementU = coinFactorizationA_->elementU();
     1301        int * indexRowU = coinFactorizationA_->indexRowU();
     1302        CoinBigIndex * startColumnU = coinFactorizationA_->startColumnU();
     1303        double slackValue = coinFactorizationA_->slackValue();
     1304        for (i=0;i<numberRowBasic;i++) {
     1305          int iRow = pivotTemp[i];
     1306          indexRowU[i]=iRow;
     1307          startColumnU[i]=i;
     1308          elementU[i]=slackValue;
     1309          numberInRow[iRow]=1;
     1310          numberInColumn[i]=1;
     1311        }
     1312        startColumnU[numberRowBasic]=numberRowBasic;
     1313        // can change for gub so redo
     1314        numberColumnBasic = numberBasic-numberRowBasic;
     1315        matrix->fillBasis(model,
     1316                          pivotTemp+numberRowBasic,
     1317                          numberColumnBasic,
     1318                          indexRowU,
     1319                          startColumnU+numberRowBasic,
     1320                          numberInRow,
     1321                          numberInColumn+numberRowBasic,
     1322                          elementU);
     1323#if 0
     1324        {
     1325          printf("%d row basic, %d column basic\n",numberRowBasic,numberColumnBasic);
     1326          for (int i=0;i<numberElements;i++)
     1327            printf("row %d col %d value %g\n",indexRowU[i],indexColumnU_[i],
     1328                   elementU[i]);
     1329        }
     1330#endif
     1331        // recompute number basic
     1332        numberBasic = numberRowBasic+numberColumnBasic;
     1333        if (numberBasic)
     1334          numberElements = startColumnU[numberBasic-1]
     1335            +numberInColumn[numberBasic-1];
     1336        else
     1337          numberElements=0;
     1338        coinFactorizationA_->setNumberElementsU(numberElements);
     1339        //saveFactorization("dump.d");
     1340        if (coinFactorizationA_->biasLU()>=3||coinFactorizationA_->numberRows()!=coinFactorizationA_->numberColumns())
     1341          coinFactorizationA_->preProcess ( 2 );
     1342        else
     1343          coinFactorizationA_->preProcess ( 3 ); // no row copy
     1344        coinFactorizationA_->factor (  );
     1345        if (coinFactorizationA_->status()==-99) {
     1346          // get more memory
     1347          coinFactorizationA_->areaFactor(2.0*coinFactorizationA_->areaFactor());
     1348        } else if (coinFactorizationA_->status()==-1&&model->numberIterations()==0&&
     1349                   coinFactorizationA_->denseThreshold()) {
     1350          // Round again without dense
     1351          coinFactorizationA_->setDenseThreshold(0);
     1352          coinFactorizationA_->setStatus(-99);
     1353        }
     1354      }
     1355      // If we get here status is 0 or -1
     1356      if (coinFactorizationA_->status() == 0) {
     1357        // We may need to tamper with order and redo - e.g. network with side
     1358        int useNumberRows = numberRows;
     1359        // **** we will also need to add test in dual steepest to do
     1360        // as we do for network
     1361        matrix->generalExpanded(model,12,useNumberRows);
     1362        const int * permuteBack = coinFactorizationA_->permuteBack();
     1363        const int * back = coinFactorizationA_->pivotColumnBack();
     1364        //int * pivotTemp = pivotColumn_.array();
     1365        //ClpDisjointCopyN ( pivotVariable, numberRows , pivotTemp  );
     1366        // Redo pivot order
     1367        for (i=0;i<numberRowBasic;i++) {
     1368          int k = pivotTemp[i];
     1369          // so rowIsBasic[k] would be permuteBack[back[i]]
     1370          pivotVariable[permuteBack[back[i]]]=k+numberColumns;
     1371        }
     1372        for (;i<useNumberRows;i++) {
     1373          int k = pivotTemp[i];
     1374          // so rowIsBasic[k] would be permuteBack[back[i]]
     1375          pivotVariable[permuteBack[back[i]]]=k;
     1376        }
     1377#if 0
     1378        if (numberSave>=0) {
     1379          numberSave=numberDense_;
     1380          memset(saveList,0,((coinFactorizationA_->numberRows()+31)>>5)*sizeof(int));
     1381          for (i=coinFactorizationA_->numberRows()-numberSave;i<coinFactorizationA_->numberRows();i++) {
     1382            int k=pivotTemp[pivotColumn_.array()[i]];
     1383            setDense(k);
     1384          }
     1385        }
     1386#endif
     1387        // Set up permutation vector
     1388        // these arrays start off as copies of permute
     1389        // (and we could use permute_ instead of pivotColumn (not back though))
     1390        ClpDisjointCopyN ( coinFactorizationA_->permute(), useNumberRows , coinFactorizationA_->pivotColumn()  );
     1391        ClpDisjointCopyN ( coinFactorizationA_->permuteBack(), useNumberRows , coinFactorizationA_->pivotColumnBack()  );
     1392#ifndef SLIM_CLP
     1393        if (networkMatrix) {
     1394          maximumPivots(CoinMax(2000,maximumPivots()));
     1395          // redo arrays
     1396          for (int iRow=0;iRow<4;iRow++) {
     1397            int length =model->numberRows()+maximumPivots();
     1398            if (iRow==3||model->objectiveAsObject()->type()>1)
     1399              length += model->numberColumns();
     1400            model->rowArray(iRow)->reserve(length);
     1401          }
     1402          // create network factorization
     1403          if (doCheck)
     1404            delete networkBasis_; // temp
     1405          networkBasis_ = new ClpNetworkBasis(model,coinFactorizationA_->numberRows(),
     1406                                              coinFactorizationA_->pivotRegion(),
     1407                                              coinFactorizationA_->permuteBack(),
     1408                                              coinFactorizationA_->startColumnU(),
     1409                                              coinFactorizationA_->numberInColumn(),
     1410                                              coinFactorizationA_->indexRowU(),
     1411                                              coinFactorizationA_->elementU());
     1412          // kill off arrays in ordinary factorization
     1413          if (!doCheck) {
     1414            coinFactorizationA_->gutsOfDestructor();
     1415            // but make sure coinFactorizationA_->numberRows() set
     1416            coinFactorizationA_->setNumberRows(model->numberRows());
     1417            coinFactorizationA_->setStatus(0);
     1418#if 0
     1419            // but put back permute arrays so odd things will work
     1420            int numberRows = model->numberRows();
     1421            pivotColumnBack_ = new int [numberRows];
     1422            permute_ = new int [numberRows];
     1423            int i;
     1424            for (i=0;i<numberRows;i++) {
     1425              pivotColumnBack_[i]=i;
     1426              permute_[i]=i;
     1427            }
     1428#endif
     1429          }
     1430        } else {
     1431#endif
     1432          // See if worth going sparse and when
     1433          coinFactorizationA_->checkSparse();
     1434#ifndef SLIM_CLP
     1435        }
     1436#endif
     1437      } else if (coinFactorizationA_->status() == -1&&(solveType==0||solveType==2)) {
     1438        // This needs redoing as it was merged coding - does not need array
     1439        int numberTotal = numberRows+numberColumns;
     1440        int * isBasic = new int [numberTotal];
     1441        int * rowIsBasic = isBasic+numberColumns;
     1442        int * columnIsBasic = isBasic;
     1443        for (i=0;i<numberTotal;i++)
     1444          isBasic[i]=-1;
     1445        for (i=0;i<numberRowBasic;i++) {
     1446          int iRow = pivotTemp[i];
     1447          rowIsBasic[iRow]=1;
     1448        }
     1449        for (;i<numberBasic;i++) {
     1450          int iColumn = pivotTemp[i];
     1451          columnIsBasic[iColumn]=1;
     1452        }
     1453        numberBasic=0;
     1454        for (i=0;i<numberRows;i++)
     1455          pivotVariable[i]=-1;
     1456        // mark as basic or non basic
     1457        const int * pivotColumn = coinFactorizationA_->pivotColumn();
     1458        for (i=0;i<numberRows;i++) {
     1459          if (rowIsBasic[i]>=0) {
     1460            if (pivotColumn[numberBasic]>=0) {
     1461              rowIsBasic[i]=pivotColumn[numberBasic];
     1462            } else {
     1463              rowIsBasic[i]=-1;
     1464              model->setRowStatus(i,ClpSimplex::superBasic);
     1465            }
     1466            numberBasic++;
     1467          }
     1468        }
     1469        for (i=0;i<numberColumns;i++) {
     1470          if (columnIsBasic[i]>=0) {
     1471            if (pivotColumn[numberBasic]>=0)
     1472              columnIsBasic[i]=pivotColumn[numberBasic];
     1473            else
     1474              columnIsBasic[i]=-1;
     1475            numberBasic++;
     1476          }
     1477        }
     1478        // leave pivotVariable in useful form for cleaning basis
     1479        int * pivotVariable = model->pivotVariable();
     1480        for (i=0;i<numberRows;i++) {
     1481          pivotVariable[i]=-1;
     1482        }
     1483
     1484        for (i=0;i<numberRows;i++) {
     1485          if (model->getRowStatus(i) == ClpSimplex::basic) {
     1486            int iPivot = rowIsBasic[i];
     1487            if (iPivot>=0)
     1488              pivotVariable[iPivot]=i+numberColumns;
     1489          }
     1490        }
     1491        for (i=0;i<numberColumns;i++) {
     1492          if (model->getColumnStatus(i) == ClpSimplex::basic) {
     1493            int iPivot = columnIsBasic[i];
     1494            if (iPivot>=0)
     1495              pivotVariable[iPivot]=i;
     1496          }
     1497        }
     1498        delete [] isBasic;
     1499        double * columnLower = model->lowerRegion();
     1500        double * columnUpper = model->upperRegion();
     1501        double * columnActivity = model->solutionRegion();
     1502        double * rowLower = model->lowerRegion(0);
     1503        double * rowUpper = model->upperRegion(0);
     1504        double * rowActivity = model->solutionRegion(0);
     1505        //redo basis - first take ALL columns out
     1506        int iColumn;
     1507        double largeValue = model->largeValue();
     1508        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1509          if (model->getColumnStatus(iColumn)==ClpSimplex::basic) {
     1510            // take out
     1511            if (!valuesPass) {
     1512              double lower=columnLower[iColumn];
     1513              double upper=columnUpper[iColumn];
     1514              double value=columnActivity[iColumn];
     1515              if (lower>-largeValue||upper<largeValue) {
     1516                if (fabs(value-lower)<fabs(value-upper)) {
     1517                  model->setColumnStatus(iColumn,ClpSimplex::atLowerBound);
     1518                  columnActivity[iColumn]=lower;
     1519                } else {
     1520                  model->setColumnStatus(iColumn,ClpSimplex::atUpperBound);
     1521                  columnActivity[iColumn]=upper;
     1522                }
     1523              } else {
     1524                model->setColumnStatus(iColumn,ClpSimplex::isFree);
     1525              }
     1526            } else {
     1527              model->setColumnStatus(iColumn,ClpSimplex::superBasic);
     1528            }
     1529          }
     1530        }
     1531        int iRow;
     1532        for (iRow=0;iRow<numberRows;iRow++) {
     1533          int iSequence=pivotVariable[iRow];
     1534          if (iSequence>=0) {
     1535            // basic
     1536            if (iSequence>=numberColumns) {
     1537              // slack in - leave
     1538              //assert (iSequence-numberColumns==iRow);
     1539            } else {
     1540              assert(model->getRowStatus(iRow)!=ClpSimplex::basic);
     1541              // put back structural
     1542              model->setColumnStatus(iSequence,ClpSimplex::basic);
     1543            }
     1544          } else {
     1545            // put in slack
     1546            model->setRowStatus(iRow,ClpSimplex::basic);
     1547          }
     1548        }
     1549        // Put back any key variables for gub
     1550        int dummy;
     1551        matrix->generalExpanded(model,1,dummy);
     1552        // signal repeat
     1553        coinFactorizationA_->setStatus(-99);
     1554        // set fixed if they are
     1555        for (iRow=0;iRow<numberRows;iRow++) {
     1556          if (model->getRowStatus(iRow)!=ClpSimplex::basic ) {
     1557            if (rowLower[iRow]==rowUpper[iRow]) {
     1558              rowActivity[iRow]=rowLower[iRow];
     1559              model->setRowStatus(iRow,ClpSimplex::isFixed);
     1560            }
     1561          }
     1562        }
     1563        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1564          if (model->getColumnStatus(iColumn)!=ClpSimplex::basic ) {
     1565            if (columnLower[iColumn]==columnUpper[iColumn]) {
     1566              columnActivity[iColumn]=columnLower[iColumn];
     1567              model->setColumnStatus(iColumn,ClpSimplex::isFixed);
     1568            }
     1569          }
     1570        }
     1571      }
     1572    }
     1573#ifndef SLIM_CLP
     1574  } else {
     1575    // network - fake factorization - do nothing
     1576    coinFactorizationA_->setStatus(0);
     1577    coinFactorizationA_->setPivots(0);
     1578  }
     1579#endif
     1580#ifndef SLIM_CLP
     1581  if (!coinFactorizationA_->status()) {
     1582    // take out part if quadratic
     1583    if (model->algorithm()==2) {
     1584      ClpObjective * obj = model->objectiveAsObject();
     1585      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
     1586      assert (quadraticObj);
     1587      CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
     1588      int numberXColumns = quadratic->getNumCols();
     1589      assert (numberXColumns<numberColumns);
     1590      int base = numberColumns-numberXColumns;
     1591      int * which = new int [numberXColumns];
     1592      int * pivotVariable = model->pivotVariable();
     1593      int * permute = pivotColumn();
     1594      int i;
     1595      int n=0;
     1596      for (i=0;i<numberRows;i++) {
     1597        int iSj = pivotVariable[i]-base;
     1598        if (iSj>=0&&iSj<numberXColumns)
     1599          which[n++]=permute[i];
     1600      }
     1601      if (n)
     1602        coinFactorizationA_->emptyRows(n,which);
     1603      delete [] which;
     1604    }
     1605  }
     1606#endif
     1607  return coinFactorizationA_->status();
     1608}
     1609/* Replaces one Column to basis,
     1610   returns 0=OK, 1=Probably OK, 2=singular, 3=no room
     1611   If checkBeforeModifying is true will do all accuracy checks
     1612   before modifying factorization.  Whether to set this depends on
     1613   speed considerations.  You could just do this on first iteration
     1614   after factorization and thereafter re-factorize
     1615   partial update already in U */
     1616int
     1617ClpFactorization::replaceColumn ( const ClpSimplex * model,
     1618                                  CoinIndexedVector * regionSparse,
     1619                                  CoinIndexedVector * tableauColumn,
     1620                                  int pivotRow,
     1621                                  double pivotCheck ,
     1622                                  bool checkBeforeModifying)
     1623{
     1624#ifndef SLIM_CLP
     1625  if (!networkBasis_) {
     1626#endif
     1627    // see if FT
     1628    if (coinFactorizationA_->forrestTomlin()) {
     1629      int returnCode= coinFactorizationA_->replaceColumn(regionSparse,
     1630                                              pivotRow,
     1631                                              pivotCheck,
     1632                                              checkBeforeModifying);
     1633      return returnCode;
     1634    } else {
     1635      return coinFactorizationA_->replaceColumnPFI(tableauColumn,
     1636                                              pivotRow,pivotCheck); // Note array
     1637    }
     1638
     1639#ifndef SLIM_CLP
     1640  } else {
     1641    if (doCheck) {
     1642      int returnCode = coinFactorizationA_->replaceColumn(regionSparse,
     1643                                                        pivotRow,
     1644                                                        pivotCheck,
     1645                                                        checkBeforeModifying);
     1646      networkBasis_->replaceColumn(regionSparse,
     1647                                   pivotRow);
     1648      return returnCode;
     1649    } else {
     1650      // increase number of pivots
     1651      coinFactorizationA_->setPivots(coinFactorizationA_->pivots()+1);
     1652      return networkBasis_->replaceColumn(regionSparse,
     1653                                   pivotRow);
     1654    }
     1655  }
     1656#endif
     1657}
     1658
     1659/* Updates one column (FTRAN) from region2
     1660   number returned is negative if no room
     1661   region1 starts as zero and is zero at end */
     1662int
     1663ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
     1664                                   CoinIndexedVector * regionSparse2)
     1665{
     1666#ifdef CLP_DEBUG
     1667  regionSparse->checkClear();
     1668#endif
     1669  if (!coinFactorizationA_->numberRows())
     1670    return 0;
     1671#ifndef SLIM_CLP
     1672  if (!networkBasis_) {
     1673#endif
     1674    coinFactorizationA_->setCollectStatistics(true);
     1675    int returnValue= coinFactorizationA_->updateColumnFT(regionSparse,
     1676                                                       regionSparse2);
     1677    coinFactorizationA_->setCollectStatistics(false);
     1678    return returnValue;
     1679#ifndef SLIM_CLP
     1680  } else {
     1681#ifdef CHECK_NETWORK
     1682    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
     1683    double * check = new double[coinFactorizationA_->numberRows()];
     1684    int returnCode = coinFactorizationA_->updateColumnFT(regionSparse,
     1685                                                       regionSparse2);
     1686    networkBasis_->updateColumn(regionSparse,save,-1);
     1687    int i;
     1688    double * array = regionSparse2->denseVector();
     1689    int * indices = regionSparse2->getIndices();
     1690    int n=regionSparse2->getNumElements();
     1691    memset(check,0,coinFactorizationA_->numberRows()*sizeof(double));
     1692    double * array2 = save->denseVector();
     1693    int * indices2 = save->getIndices();
     1694    int n2=save->getNumElements();
     1695    assert (n==n2);
     1696    if (save->packedMode()) {
     1697      for (i=0;i<n;i++) {
     1698        check[indices[i]]=array[i];
     1699      }
     1700      for (i=0;i<n;i++) {
     1701        double value2 = array2[i];
     1702        assert (check[indices2[i]]==value2);
     1703      }
     1704    } else {
     1705      int numberRows = coinFactorizationA_->numberRows();
     1706      for (i=0;i<numberRows;i++) {
     1707        double value1 = array[i];
     1708        double value2 = array2[i];
     1709        assert (value1==value2);
     1710      }
     1711    }
     1712    delete save;
     1713    delete [] check;
     1714    return returnCode;
     1715#else
     1716    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
     1717    return 1;
     1718#endif
     1719  }
     1720#endif
     1721}
     1722/* Updates one column (FTRAN) from region2
     1723   number returned is negative if no room
     1724   region1 starts as zero and is zero at end */
     1725int
     1726ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
     1727                                 CoinIndexedVector * regionSparse2,
     1728                                 bool noPermute) const
     1729{
     1730#ifdef CLP_DEBUG
     1731  if (!noPermute)
     1732    regionSparse->checkClear();
     1733#endif
     1734  if (!coinFactorizationA_->numberRows())
     1735    return 0;
     1736#ifndef SLIM_CLP
     1737  if (!networkBasis_) {
     1738#endif
     1739    coinFactorizationA_->setCollectStatistics(true);
     1740    int returnValue= coinFactorizationA_->updateColumn(regionSparse,
     1741                                                     regionSparse2,
     1742                                                     noPermute);
     1743    coinFactorizationA_->setCollectStatistics(false);
     1744    return returnValue;
     1745#ifndef SLIM_CLP
     1746  } else {
     1747#ifdef CHECK_NETWORK
     1748    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
     1749    double * check = new double[coinFactorizationA_->numberRows()];
     1750    int returnCode = coinFactorizationA_->updateColumn(regionSparse,
     1751                                                     regionSparse2,
     1752                                                     noPermute);
     1753    networkBasis_->updateColumn(regionSparse,save,-1);
     1754    int i;
     1755    double * array = regionSparse2->denseVector();
     1756    int * indices = regionSparse2->getIndices();
     1757    int n=regionSparse2->getNumElements();
     1758    memset(check,0,coinFactorizationA_->numberRows()*sizeof(double));
     1759    double * array2 = save->denseVector();
     1760    int * indices2 = save->getIndices();
     1761    int n2=save->getNumElements();
     1762    assert (n==n2);
     1763    if (save->packedMode()) {
     1764      for (i=0;i<n;i++) {
     1765        check[indices[i]]=array[i];
     1766      }
     1767      for (i=0;i<n;i++) {
     1768        double value2 = array2[i];
     1769        assert (check[indices2[i]]==value2);
     1770      }
     1771    } else {
     1772      int numberRows = coinFactorizationA_->numberRows();
     1773      for (i=0;i<numberRows;i++) {
     1774        double value1 = array[i];
     1775        double value2 = array2[i];
     1776        assert (value1==value2);
     1777      }
     1778    }
     1779    delete save;
     1780    delete [] check;
     1781    return returnCode;
     1782#else
     1783    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
     1784    return 1;
     1785#endif
     1786  }
     1787#endif
     1788}
     1789/* Updates one column (FTRAN) from region2
     1790   number returned is negative if no room
     1791   region1 starts as zero and is zero at end */
     1792int
     1793ClpFactorization::updateColumnForDebug ( CoinIndexedVector * regionSparse,
     1794                                 CoinIndexedVector * regionSparse2,
     1795                                 bool noPermute) const
     1796{
     1797  if (!noPermute)
     1798    regionSparse->checkClear();
     1799  if (!coinFactorizationA_->numberRows())
     1800    return 0;
     1801  coinFactorizationA_->setCollectStatistics(false);
     1802  int returnValue= coinFactorizationA_->updateColumn(regionSparse,
     1803                                                   regionSparse2,
     1804                                                   noPermute);
     1805  return returnValue;
     1806}
     1807/* Updates one column (BTRAN) from region2
     1808   region1 starts as zero and is zero at end */
     1809int
     1810ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
     1811                          CoinIndexedVector * regionSparse2) const
     1812{
     1813  if (!coinFactorizationA_->numberRows())
     1814    return 0;
     1815#ifndef SLIM_CLP
     1816  if (!networkBasis_) {
     1817#endif
     1818    coinFactorizationA_->setCollectStatistics(true);
     1819    return coinFactorizationA_->updateColumnTranspose(regionSparse,
     1820                                                    regionSparse2);
     1821    coinFactorizationA_->setCollectStatistics(false);
     1822#ifndef SLIM_CLP
     1823  } else {
     1824#ifdef CHECK_NETWORK
     1825    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
     1826    double * check = new double[coinFactorizationA_->numberRows()];
     1827    int returnCode = coinFactorizationA_->updateColumnTranspose(regionSparse,
     1828                                                              regionSparse2);
     1829    networkBasis_->updateColumnTranspose(regionSparse,save);
     1830    int i;
     1831    double * array = regionSparse2->denseVector();
     1832    int * indices = regionSparse2->getIndices();
     1833    int n=regionSparse2->getNumElements();
     1834    memset(check,0,coinFactorizationA_->numberRows()*sizeof(double));
     1835    double * array2 = save->denseVector();
     1836    int * indices2 = save->getIndices();
     1837    int n2=save->getNumElements();
     1838    assert (n==n2);
     1839    if (save->packedMode()) {
     1840      for (i=0;i<n;i++) {
     1841        check[indices[i]]=array[i];
     1842      }
     1843      for (i=0;i<n;i++) {
     1844        double value2 = array2[i];
     1845        assert (check[indices2[i]]==value2);
     1846      }
     1847    } else {
     1848      int numberRows = coinFactorizationA_->numberRows();
     1849      for (i=0;i<numberRows;i++) {
     1850        double value1 = array[i];
     1851        double value2 = array2[i];
     1852        assert (value1==value2);
     1853      }
     1854    }
     1855    delete save;
     1856    delete [] check;
     1857    return returnCode;
     1858#else
     1859    return networkBasis_->updateColumnTranspose(regionSparse,regionSparse2);
     1860#endif
     1861  }
     1862#endif
     1863}
     1864/* makes a row copy of L for speed and to allow very sparse problems */
     1865void
     1866ClpFactorization::goSparse()
     1867{
     1868#ifndef SLIM_CLP
     1869  if (!networkBasis_)
     1870#endif
     1871    coinFactorizationA_->goSparse();
     1872}
     1873// Cleans up i.e. gets rid of network basis
     1874void
     1875ClpFactorization::cleanUp()
     1876{
     1877#ifndef SLIM_CLP
     1878  delete networkBasis_;
     1879  networkBasis_=NULL;
     1880#endif
     1881  coinFactorizationA_->resetStatistics();
     1882}
     1883/// Says whether to redo pivot order
     1884bool
     1885ClpFactorization::needToReorder() const
     1886{
     1887#ifdef CHECK_NETWORK
     1888  return true;
     1889#endif
     1890#ifndef SLIM_CLP
     1891  if (!networkBasis_)
     1892#endif
     1893    return true;
     1894#ifndef SLIM_CLP
     1895  else
     1896    return false;
     1897#endif
     1898}
     1899// Get weighted row list
     1900void
     1901ClpFactorization::getWeights(int * weights) const
     1902{
     1903#ifndef SLIM_CLP
     1904  if (networkBasis_) {
     1905    // Network - just unit
     1906    int numberRows = coinFactorizationA_->numberRows();
     1907    for (int i=0;i<numberRows;i++)
     1908      weights[i]=1;
     1909    return;
     1910  }
     1911#endif
     1912  int * numberInRow = coinFactorizationA_->numberInRow();
     1913  int * numberInColumn = coinFactorizationA_->numberInColumn();
     1914  int * permuteBack = coinFactorizationA_->pivotColumnBack();
     1915  int * indexRowU = coinFactorizationA_->indexRowU();
     1916  const CoinBigIndex * startColumnU = coinFactorizationA_->startColumnU();
     1917  const CoinBigIndex * startRowL = coinFactorizationA_->startRowL();
     1918  int numberRows = coinFactorizationA_->numberRows();
     1919  if (!startRowL||!coinFactorizationA_->numberInRow()) {
     1920    int * temp = new int[numberRows];
     1921    memset(temp,0,numberRows*sizeof(int));
     1922    int i;
     1923    for (i=0;i<numberRows;i++) {
     1924      // one for pivot
     1925      temp[i]++;
     1926      CoinBigIndex j;
     1927      for (j=startColumnU[i];j<startColumnU[i]+numberInColumn[i];j++) {
     1928        int iRow=indexRowU[j];
     1929        temp[iRow]++;
     1930      }
     1931    }
     1932    CoinBigIndex * startColumnL = coinFactorizationA_->startColumnL();
     1933    int * indexRowL = coinFactorizationA_->indexRowL();
     1934    int numberL = coinFactorizationA_->numberL();
     1935    CoinBigIndex baseL = coinFactorizationA_->baseL();
     1936    for (i=baseL;i<baseL+numberL;i++) {
     1937      CoinBigIndex j;
     1938      for (j=startColumnL[i];j<startColumnL[i+1];j++) {
     1939        int iRow = indexRowL[j];
     1940        temp[iRow]++;
     1941      }
     1942    }
     1943    for (i=0;i<numberRows;i++) {
     1944      int number = temp[i];
     1945      int iPermute = permuteBack[i];
     1946      weights[iPermute]=number;
     1947    }
     1948    delete [] temp;
     1949  } else {
     1950    int i;
     1951    for (i=0;i<numberRows;i++) {
     1952      int number = startRowL[i+1]-startRowL[i]+numberInRow[i]+1;
     1953      //number = startRowL[i+1]-startRowL[i]+1;
     1954      //number = numberInRow[i]+1;
     1955      int iPermute = permuteBack[i];
     1956      weights[iPermute]=number;
     1957    }
     1958  }
     1959}
     1960#endif
  • trunk/Clp/src/ClpFactorization.hpp

    r1055 r1147  
    1616    a genuine ClpNetworkBasis object
    1717*/
    18 
    19 class ClpFactorization : public CoinFactorization {
     18class ClpFactorization
     19#ifndef CLP_MULTIPLE_FACTORIZATIONS   
     20  : public CoinFactorization
     21#endif
     22{
     23
     24  //friend class CoinFactorization;
    2025 
    2126public:
     
    9398                              CoinIndexedVector * regionSparse2) const;
    9499  //@}
     100#ifdef CLP_MULTIPLE_FACTORIZATIONS   
     101  /**@name Lifted from CoinFactorization */
     102  //@{
     103  /// Total number of elements in factorization
     104  inline int numberElements (  ) const {
     105    if (coinFactorizationA_) return coinFactorizationA_->numberElements(); else return 0 ;
     106  }
     107  /// Returns address of permute region
     108  inline int *permute (  ) const {
     109    if (coinFactorizationA_) return coinFactorizationA_->permute(); else return NULL ;
     110  }
     111  /// Returns address of pivotColumn region (also used for permuting)
     112  inline int *pivotColumn (  ) const {
     113    if (coinFactorizationA_) return coinFactorizationA_->pivotColumn(); else return NULL ;
     114  }
     115  /// Maximum number of pivots between factorizations
     116  inline int maximumPivots (  ) const {
     117    if (coinFactorizationA_) return coinFactorizationA_->maximumPivots();  else return 0 ;
     118  }
     119  /// Set maximum number of pivots between factorizations
     120  inline void maximumPivots (  int value) {
     121    if (coinFactorizationA_) coinFactorizationA_->maximumPivots(value);
     122  }
     123  /// Returns number of pivots since factorization
     124  inline int pivots (  ) const {
     125    if (coinFactorizationA_) return coinFactorizationA_->pivots(); else return 0 ;
     126  }
     127  /// Whether larger areas needed
     128  inline double areaFactor (  ) const {
     129    if (coinFactorizationA_) return coinFactorizationA_->areaFactor(); else return 0.0 ;
     130  }
     131  /// Set whether larger areas needed
     132  inline void areaFactor ( double value) {
     133    if (coinFactorizationA_) coinFactorizationA_->areaFactor(value);
     134  }
     135  /// Zero tolerance
     136  inline double zeroTolerance (  ) const {
     137    if (coinFactorizationA_) return coinFactorizationA_->zeroTolerance();  else return 1.0e-13 ;
     138  }
     139  /// Set zero tolerance
     140  inline void zeroTolerance (  double value) {
     141    if (coinFactorizationA_) coinFactorizationA_->zeroTolerance(value);
     142  }
     143  /**  get sparse threshold */
     144  inline int sparseThreshold ( ) const
     145  { if (coinFactorizationA_) return coinFactorizationA_->sparseThreshold(); else return 0 ;}
     146  /**  Set sparse threshold */
     147  inline void sparseThreshold ( int value)
     148  { if (coinFactorizationA_) coinFactorizationA_->sparseThreshold(value); }
     149  /// Returns status
     150  inline int status (  ) const {
     151    if (coinFactorizationA_) return coinFactorizationA_->status(); else return 0 ;
     152  }
     153  /// Returns number of dense rows
     154  inline int numberDense() const
     155  { if (coinFactorizationA_) return coinFactorizationA_->numberDense(); else return 0 ;}
     156#if 0
     157  /// Returns number in U area
     158  inline CoinBigIndex numberElementsU (  ) const {
     159    if (coinFactorizationA_) return coinFactorizationA_->numberElementsU(); else return 0 ;
     160  }
     161  /// Returns number in L area
     162  inline CoinBigIndex numberElementsL (  ) const {
     163    if (coinFactorizationA_) return coinFactorizationA_->numberElementsL(); else return 0 ;
     164  }
     165  /// Returns number in R area
     166  inline CoinBigIndex numberElementsR (  ) const {
     167    if (coinFactorizationA_) return coinFactorizationA_->numberElementsR(); else return 0 ;
     168  }
     169#endif
     170  inline bool timeToRefactorize() const
     171  {
     172    if (coinFactorizationA_) {
     173      return (coinFactorizationA_->pivots()*3>coinFactorizationA_->maximumPivots()*2&&
     174              coinFactorizationA_->numberElementsR()*3>(coinFactorizationA_->numberElementsL()+
     175                                                        coinFactorizationA_->numberElementsU())*2+1000&&
     176              !coinFactorizationA_->numberDense());
     177    } else {
     178      return false;
     179    }
     180  }
     181  /// Level of detail of messages
     182  inline int messageLevel (  ) const {
     183    if (coinFactorizationA_) return coinFactorizationA_->messageLevel();  else return 0 ;
     184  }
     185  /// Set level of detail of messages
     186  inline void messageLevel (  int value) {
     187    if (coinFactorizationA_) coinFactorizationA_->messageLevel(value);
     188  }
     189  /// Get rid of all memory
     190  inline void clearArrays()
     191  { if (coinFactorizationA_) coinFactorizationA_->clearArrays();}
     192  /// Number of Rows after factorization
     193  inline int numberRows (  ) const {
     194    if (coinFactorizationA_) return coinFactorizationA_->numberRows(); else return 0 ;
     195  }
     196  /// Gets dense threshold
     197  inline int denseThreshold() const
     198  { if (coinFactorizationA_) return coinFactorizationA_->denseThreshold(); else return 0 ;}
     199  /// Sets dense threshold
     200  inline void setDenseThreshold(int value)
     201  { if (coinFactorizationA_) coinFactorizationA_->setDenseThreshold(value);}
     202  /// Pivot tolerance
     203  inline double pivotTolerance (  ) const {
     204    if (coinFactorizationA_) return coinFactorizationA_->pivotTolerance();  else return 1.0e-8 ;
     205  }
     206  /// Set pivot tolerance
     207  inline void pivotTolerance (  double value) {
     208    if (coinFactorizationA_) coinFactorizationA_->pivotTolerance(value);
     209  }
     210  /// Allows change of pivot accuracy check 1.0 == none >1.0 relaxed
     211  inline void relaxAccuracyCheck(double value)
     212  { if (coinFactorizationA_) coinFactorizationA_->relaxAccuracyCheck(value);}
     213  /** Array persistence flag
     214      If 0 then as now (delete/new)
     215      1 then only do arrays if bigger needed
     216      2 as 1 but give a bit extra if bigger needed
     217  */
     218  inline int persistenceFlag() const
     219  { if (coinFactorizationA_) return coinFactorizationA_->persistenceFlag(); else return 0 ;}
     220  inline void setPersistenceFlag(int value)
     221  { if (coinFactorizationA_) coinFactorizationA_->setPersistenceFlag(value);}
     222  /// Delete all stuff (leaves as after CoinFactorization())
     223  inline void almostDestructor()
     224  { if (coinFactorizationA_) coinFactorizationA_->almostDestructor(); }
     225  /// Returns areaFactor but adjusted for dense
     226  inline double adjustedAreaFactor() const
     227  { if (coinFactorizationA_) return coinFactorizationA_->adjustedAreaFactor(); else return 0.0 ; }
     228  inline void setBiasLU(int value)
     229  { if (coinFactorizationA_) coinFactorizationA_->setBiasLU(value);}
     230  /// true if Forrest Tomlin update, false if PFI
     231  inline void setForrestTomlin(bool value)
     232  { if (coinFactorizationA_) coinFactorizationA_->setForrestTomlin(value);}
     233  /// Sets default values
     234  inline void setDefaultValues() {
     235    if (coinFactorizationA_) {
     236      coinFactorizationA_->increasingRows(2);
     237      // row activities have negative sign
     238      coinFactorizationA_->slackValue(-1.0);
     239      coinFactorizationA_->zeroTolerance(1.0e-13);
     240    }
     241  }
     242#else
     243  inline bool timeToRefactorize() const
     244  {
     245    return (pivots()*3>maximumPivots()*2&&
     246            numberElementsR()*3>(numberElementsL()+numberElementsU())*2+1000&&
     247            !numberDense());
     248  }
     249  /// Sets default values
     250  inline void setDefaultValues() {
     251    increasingRows(2);
     252    // row activities have negative sign
     253    slackValue(-1.0);
     254    zeroTolerance(1.0e-13);
     255  }
     256#endif
     257  //@}
    95258   
    96259  /**@name other stuff */
     
    124287  ClpNetworkBasis * networkBasis_;
    125288#endif
     289#ifdef CLP_MULTIPLE_FACTORIZATIONS   
     290  /// Pointer to CoinFactorization
     291  CoinFactorization * coinFactorizationA_;
     292#endif
    126293  //@}
    127294};
  • trunk/Clp/src/ClpMain.cpp

    r1111 r1147  
    1111#include <iostream>
    1212
     13int boundary_sort=1000;
     14int boundary_sort2=1000;
     15int boundary_sort3=10000;
    1316
    1417#include "CoinPragma.hpp"
  • trunk/Clp/src/ClpModel.cpp

    r1141 r1147  
    8080  defaultHandler_(true),
    8181  rowNames_(),
    82   columnNames_()
     82  columnNames_(),
    8383#else
    84   defaultHandler_(true)
    85 #endif
     84  defaultHandler_(true),
     85#endif
     86  maximumColumns_(-1),
     87  maximumRows_(-1),
     88  savedRowScale_(NULL),
     89  savedColumnScale_(NULL)
    8690{
    8791  intParam_[ClpMaxNumIteration] = 2147483647;
     
    118122    handler_ = NULL;
    119123  }
    120   gutsOfDelete();
    121 }
    122 void ClpModel::gutsOfDelete()
    123 {
    124   delete [] rowActivity_;
    125   rowActivity_=NULL;
    126   delete [] columnActivity_;
    127   columnActivity_=NULL;
    128   delete [] dual_;
    129   dual_=NULL;
    130   delete [] reducedCost_;
    131   reducedCost_=NULL;
    132   delete [] rowLower_;
    133   delete [] rowUpper_;
    134   delete [] rowObjective_;
    135   rowLower_=NULL;
    136   rowUpper_=NULL;
    137   rowObjective_=NULL;
    138   delete [] columnLower_;
    139   delete [] columnUpper_;
    140   delete objective_;
    141   columnLower_=NULL;
    142   columnUpper_=NULL;
    143   objective_=NULL;
     124  gutsOfDelete(0);
     125}
     126// Does most of deletion (0 = all, 1 = most)
     127void
     128ClpModel::gutsOfDelete(int type)
     129{
     130  if (!type||!permanentArrays()) {
     131    maximumRows_=-1;
     132    maximumColumns_ = -1;
     133    delete [] rowActivity_;
     134    rowActivity_=NULL;
     135    delete [] columnActivity_;
     136    columnActivity_=NULL;
     137    delete [] dual_;
     138    dual_=NULL;
     139    delete [] reducedCost_;
     140    reducedCost_=NULL;
     141    delete [] rowLower_;
     142    delete [] rowUpper_;
     143    delete [] rowObjective_;
     144    rowLower_=NULL;
     145    rowUpper_=NULL;
     146    rowObjective_=NULL;
     147    delete [] columnLower_;
     148    delete [] columnUpper_;
     149    delete objective_;
     150    columnLower_=NULL;
     151    columnUpper_=NULL;
     152    objective_=NULL;
     153    delete [] savedRowScale_;
     154    if (rowScale_==savedRowScale_)
     155      rowScale_=NULL;
     156    savedRowScale_ = NULL;
     157    delete [] savedColumnScale_;
     158    if (columnScale_==savedColumnScale_)
     159      columnScale_=NULL;
     160    savedColumnScale_ = NULL;
     161    delete [] rowScale_;
     162    rowScale_ = NULL;
     163    delete [] columnScale_;
     164    columnScale_ = NULL;
     165    delete [] integerType_;
     166    integerType_ = NULL;
     167    delete [] status_;
     168    status_=NULL;
     169    delete eventHandler_;
     170    eventHandler_=NULL;
     171  }
     172  whatsChanged_=0;
    144173  delete matrix_;
    145174  matrix_=NULL;
     
    148177  delete [] ray_;
    149178  ray_ = NULL;
    150   delete [] rowScale_;
    151   rowScale_ = NULL;
    152   delete [] columnScale_;
    153   columnScale_ = NULL;
    154   delete [] integerType_;
    155   integerType_ = NULL;
    156   delete [] status_;
    157   status_=NULL;
    158   delete eventHandler_;
    159   eventHandler_=NULL;
    160   whatsChanged_=0;
    161179  specialOptions_ = 0;
     180}
     181void
     182ClpModel::setRowScale(double * scale)
     183{
     184  if (!savedRowScale_) {
     185    delete [] (double *) rowScale_;
     186    rowScale_ = scale;
     187  } else {
     188    assert (!scale);
     189    rowScale_ = NULL;
     190  }
     191}
     192void
     193ClpModel::setColumnScale(double * scale)
     194{
     195  if (!savedColumnScale_) {
     196    delete [] (double *) columnScale_;
     197    columnScale_ = scale;
     198  } else {
     199    assert (!scale);
     200    columnScale_ = NULL;
     201  }
    162202}
    163203//#############################################################################
     
    187227  // Save specialOptions
    188228  int saveOptions = specialOptions_;
    189   gutsOfDelete();
     229  gutsOfDelete(0);
    190230  specialOptions_ = saveOptions;
    191231  eventHandler_ = handler;
     
    597637  optimizationDirection_(rhs.optimizationDirection_),
    598638  numberRows_(rhs.numberRows_),
    599   numberColumns_(rhs.numberColumns_)
     639  numberColumns_(rhs.numberColumns_),
     640  specialOptions_(rhs.specialOptions_),
     641  maximumColumns_(-1),
     642  maximumRows_(-1),
     643  savedRowScale_(NULL),
     644  savedColumnScale_(NULL)
    600645{
    601646  gutsOfCopy(rhs);
     
    604649    // really do scaling
    605650    scalingFlag_=scalingMode;
    606     delete [] rowScale_;
    607     rowScale_ = NULL;
    608     delete [] columnScale_;
    609     columnScale_ = NULL;
     651    setRowScale(NULL);
     652    setColumnScale(NULL);
    610653    delete rowCopy_; // in case odd
    611654    rowCopy_=NULL;
     
    628671  if (this != &rhs) {
    629672    if (defaultHandler_) {
    630       delete handler_;
    631       handler_ = NULL;
    632     }
    633     gutsOfDelete();
     673      //delete handler_;
     674      //handler_ = NULL;
     675    }
     676    gutsOfDelete(1);
    634677    optimizationDirection_ = rhs.optimizationDirection_;
    635678    numberRows_ = rhs.numberRows_;
    636679    numberColumns_ = rhs.numberColumns_;
    637     gutsOfCopy(rhs);
     680    gutsOfCopy(rhs,-1);
    638681  }
    639682  return *this;
    640683}
    641 // Does most of copying
    642 void
    643 ClpModel::gutsOfCopy(const ClpModel & rhs, bool trueCopy)
     684/* Does most of copying
     685   If trueCopy 0 then just points to arrays
     686   If -1 leaves as much as possible */
     687void
     688ClpModel::gutsOfCopy(const ClpModel & rhs, int trueCopy)
    644689{
    645690  defaultHandler_ = rhs.defaultHandler_;
    646   if (defaultHandler_)
    647     handler_ = new CoinMessageHandler(*rhs.handler_);
    648    else
    649     handler_ = rhs.handler_;
    650   eventHandler_ = rhs.eventHandler_->clone();
    651   randomNumberGenerator_ = rhs.randomNumberGenerator_;
    652   messages_ = rhs.messages_;
    653   coinMessages_ = rhs.coinMessages_;
     691  if (trueCopy>=0) {
     692    if (defaultHandler_)
     693      handler_ = new CoinMessageHandler(*rhs.handler_);
     694    else
     695      handler_ = rhs.handler_;
     696    eventHandler_ = rhs.eventHandler_->clone();
     697    randomNumberGenerator_ = rhs.randomNumberGenerator_;
     698    messages_ = rhs.messages_;
     699    coinMessages_ = rhs.coinMessages_;
     700  } else {
     701    if (!eventHandler_&&rhs.eventHandler_)
     702      eventHandler_ = rhs.eventHandler_->clone();
     703  }
    654704  intParam_[ClpMaxNumIteration] = rhs.intParam_[ClpMaxNumIteration];
    655705  intParam_[ClpMaxNumIterationHotStart] =
     
    693743#endif
    694744    numberThreads_ = rhs.numberThreads_;
    695     integerType_ = CoinCopyOfArray(rhs.integerType_,numberColumns_);
    696     if (rhs.rowActivity_) {
    697       rowActivity_=new double[numberRows_];
    698       columnActivity_=new double[numberColumns_];
    699       dual_=new double[numberRows_];
    700       reducedCost_=new double[numberColumns_];
    701       ClpDisjointCopyN ( rhs.rowActivity_, numberRows_ ,
    702                           rowActivity_);
    703       ClpDisjointCopyN ( rhs.columnActivity_, numberColumns_ ,
    704                           columnActivity_);
    705       ClpDisjointCopyN ( rhs.dual_, numberRows_ ,
    706                           dual_);
    707       ClpDisjointCopyN ( rhs.reducedCost_, numberColumns_ ,
    708                           reducedCost_);
     745    if (maximumRows_<0) {
     746      specialOptions_ &= ~65536;
     747      savedRowScale_ = NULL;
     748      savedColumnScale_ = NULL;
     749      integerType_ = CoinCopyOfArray(rhs.integerType_,numberColumns_);
     750      if (rhs.rowActivity_) {
     751        rowActivity_=new double[numberRows_];
     752        columnActivity_=new double[numberColumns_];
     753        dual_=new double[numberRows_];
     754        reducedCost_=new double[numberColumns_];
     755        ClpDisjointCopyN ( rhs.rowActivity_, numberRows_ ,
     756                           rowActivity_);
     757        ClpDisjointCopyN ( rhs.columnActivity_, numberColumns_ ,
     758                           columnActivity_);
     759        ClpDisjointCopyN ( rhs.dual_, numberRows_ ,
     760                           dual_);
     761        ClpDisjointCopyN ( rhs.reducedCost_, numberColumns_ ,
     762                           reducedCost_);
     763      } else {
     764        rowActivity_=NULL;
     765        columnActivity_=NULL;
     766        dual_=NULL;
     767        reducedCost_=NULL;
     768      }
     769      rowLower_ = ClpCopyOfArray ( rhs.rowLower_, numberRows_ );
     770      rowUpper_ = ClpCopyOfArray ( rhs.rowUpper_, numberRows_ );
     771      columnLower_ = ClpCopyOfArray ( rhs.columnLower_, numberColumns_ );
     772      int scaleLength = ((specialOptions_&131072)==0) ? 1 : 2;
     773      columnUpper_ = ClpCopyOfArray ( rhs.columnUpper_, numberColumns_ );
     774      rowScale_ = ClpCopyOfArray(rhs.rowScale_,numberRows_*scaleLength);
     775      columnScale_ = ClpCopyOfArray(rhs.columnScale_,numberColumns_*scaleLength);
     776      if (rhs.objective_)
     777        objective_  = rhs.objective_->clone();
     778      else
     779        objective_ = NULL;
     780      rowObjective_ = ClpCopyOfArray ( rhs.rowObjective_, numberRows_ );
     781      status_ = ClpCopyOfArray( rhs.status_,numberColumns_+numberRows_);
     782      ray_ = NULL;
     783      if (problemStatus_==1&&!secondaryStatus_)
     784        ray_ = ClpCopyOfArray (rhs.ray_,numberRows_);
     785      else if (problemStatus_==2)
     786        ray_ = ClpCopyOfArray (rhs.ray_,numberColumns_);
     787      if (rhs.rowCopy_) {
     788        rowCopy_ = rhs.rowCopy_->clone();
     789      } else {
     790        rowCopy_=NULL;
     791      }
     792      matrix_=NULL;
     793      if (rhs.matrix_) {
     794        matrix_ = rhs.matrix_->clone();
     795      }
    709796    } else {
    710       rowActivity_=NULL;
    711       columnActivity_=NULL;
    712       dual_=NULL;
    713       reducedCost_=NULL;
    714     }
    715     rowLower_ = ClpCopyOfArray ( rhs.rowLower_, numberRows_ );
    716     rowUpper_ = ClpCopyOfArray ( rhs.rowUpper_, numberRows_ );
    717     columnLower_ = ClpCopyOfArray ( rhs.columnLower_, numberColumns_ );
    718     int scaleLength = ((specialOptions_&131072)==0) ? 1 : 2;
    719     columnUpper_ = ClpCopyOfArray ( rhs.columnUpper_, numberColumns_ );
    720     rowScale_ = ClpCopyOfArray(rhs.rowScale_,numberRows_*scaleLength);
    721     columnScale_ = ClpCopyOfArray(rhs.columnScale_,numberColumns_*scaleLength);
    722     if (rhs.objective_)
    723       objective_  = rhs.objective_->clone();
    724     else
    725       objective_ = NULL;
    726     rowObjective_ = ClpCopyOfArray ( rhs.rowObjective_, numberRows_ );
    727     status_ = ClpCopyOfArray( rhs.status_,numberColumns_+numberRows_);
    728     ray_ = NULL;
    729     if (problemStatus_==1&&!secondaryStatus_)
    730       ray_ = ClpCopyOfArray (rhs.ray_,numberRows_);
    731     else if (problemStatus_==2)
    732       ray_ = ClpCopyOfArray (rhs.ray_,numberColumns_);
    733     if (rhs.rowCopy_) {
    734       rowCopy_ = rhs.rowCopy_->clone();
    735     } else {
    736       rowCopy_=NULL;
    737     }
    738     matrix_=NULL;
    739     if (rhs.matrix_) {
    740       matrix_ = rhs.matrix_->clone();
     797      // This already has arrays - just copy
     798      startPermanentArrays();
     799      if (rhs.integerType_) {
     800        assert (integerType_);
     801        ClpDisjointCopyN(rhs.integerType_,numberColumns_,integerType_);
     802      } else {
     803        integerType_ = NULL;
     804      }
     805      if (rhs.rowActivity_) {
     806        ClpDisjointCopyN ( rhs.rowActivity_, numberRows_ ,
     807                           rowActivity_);
     808        ClpDisjointCopyN ( rhs.columnActivity_, numberColumns_ ,
     809                           columnActivity_);
     810        ClpDisjointCopyN ( rhs.dual_, numberRows_ ,
     811                           dual_);
     812        ClpDisjointCopyN ( rhs.reducedCost_, numberColumns_ ,
     813                           reducedCost_);
     814      } else {
     815        rowActivity_=NULL;
     816        columnActivity_=NULL;
     817        dual_=NULL;
     818        reducedCost_=NULL;
     819      }
     820      ClpDisjointCopyN ( rhs.rowLower_, numberRows_,rowLower_ );
     821      ClpDisjointCopyN ( rhs.rowUpper_, numberRows_,rowUpper_ );
     822      ClpDisjointCopyN ( rhs.columnLower_, numberColumns_,columnLower_ );
     823      assert ((specialOptions_&131072)==0);
     824      ClpDisjointCopyN ( rhs.columnUpper_, numberColumns_,columnUpper_ );
     825      if (rhs.objective_) {
     826        abort(); //check if same
     827        objective_  = rhs.objective_->clone();
     828      } else {
     829        objective_ = NULL;
     830      }
     831      assert (!rhs.rowObjective_);
     832      ClpDisjointCopyN( rhs.status_,numberColumns_+numberRows_,status_);
     833      ray_ = NULL;
     834      if (problemStatus_==1&&!secondaryStatus_)
     835        ray_ = ClpCopyOfArray (rhs.ray_,numberRows_);
     836      else if (problemStatus_==2)
     837        ray_ = ClpCopyOfArray (rhs.ray_,numberColumns_);
     838      assert (!ray_);
     839      delete rowCopy_;
     840      if (rhs.rowCopy_) {
     841        rowCopy_ = rhs.rowCopy_->clone();
     842      } else {
     843        rowCopy_=NULL;
     844      }
     845      delete matrix_;
     846      matrix_=NULL;
     847      if (rhs.matrix_) {
     848        matrix_ = rhs.matrix_->clone();
     849      }
     850      if (rhs.savedRowScale_) {
     851        assert (savedRowScale_);
     852        assert (!rowScale_);
     853        ClpDisjointCopyN(rhs.savedRowScale_,2*maximumRows_,savedRowScale_);
     854        ClpDisjointCopyN(rhs.savedColumnScale_,2*maximumColumns_,savedColumnScale_);
     855      } else {
     856        assert (!savedRowScale_);
     857        if (rowScale_) {
     858          ClpDisjointCopyN(rhs.rowScale_,numberRows_,rowScale_);
     859          ClpDisjointCopyN(rhs.columnScale_,numberColumns_,columnScale_);
     860        } else {
     861          rowScale_ =NULL;
     862          columnScale_=NULL;
     863        }
     864      }
     865      abort(); // look at resizeDouble and resize
    741866    }
    742867  } else {
     868    savedRowScale_ = rhs.savedRowScale_;
     869    savedColumnScale_ = rhs.savedColumnScale_;
    743870    rowActivity_ = rhs.rowActivity_;
    744871    columnActivity_ = rhs.columnActivity_;
     
    766893  }
    767894}
     895// Copy contents - resizing if necessary - otherwise re-use memory
     896void
     897ClpModel::copy(const ClpMatrixBase * from, ClpMatrixBase * & to)
     898{
     899  assert (from);
     900  assert (to);
     901  const ClpPackedMatrix * matrixFrom = (dynamic_cast<const ClpPackedMatrix*>(from));
     902  ClpPackedMatrix * matrixTo = (dynamic_cast< ClpPackedMatrix*>(to));
     903  if (matrixFrom&&matrixTo) {
     904    matrixTo->copy(matrixFrom);
     905  } else {
     906    delete to;
     907    to =  from->clone();
     908  }
     909#if 0
     910    delete modelPtr_->matrix_;
     911    if (continuousModel_->matrix_) {
     912      modelPtr_->matrix_ = continuousModel_->matrix_->clone();
     913    } else {
     914      modelPtr_->matrix_=NULL;
     915    }
     916#endif
     917}
    768918/* Borrow model.  This is so we dont have to copy large amounts
    769919   of data around.  It assumes a derived class wants to overwrite
     
    776926    handler_ = NULL;
    777927  }
    778   gutsOfDelete();
     928  gutsOfDelete(1);
    779929  optimizationDirection_ = rhs.optimizationDirection_;
    780930  numberRows_ = rhs.numberRows_;
     
    782932  delete [] rhs.ray_;
    783933  rhs.ray_=NULL;
    784   gutsOfCopy(rhs,false);
     934  gutsOfCopy(rhs,0);
     935  specialOptions_ = rhs.specialOptions_ & ~65536;
     936  savedRowScale_=NULL;
     937  savedColumnScale_=NULL;
    785938}
    786939// Return model - nulls all arrays so can be deleted safely
     
    807960  otherModel.ray_ = ray_;
    808961  ray_ = NULL;
     962  rowScale_ = NULL;
     963  columnScale_ = NULL;
    809964  //rowScale_=NULL;
    810965  //columnScale_=NULL;
     
    9151070                      bool createArray)
    9161071{
    917   if ((array||createArray)&&size!=newSize) {
     1072  if ((array||createArray)&&size<newSize) {
    9181073    int i;
    9191074    double * newArray = new double[newSize];
     
    10361191  reducedCost_ = resizeDouble(reducedCost_,numberColumns_,
    10371192                              newNumberColumns,0.0,true);
     1193  savedRowScale_ = resizeDouble(savedRowScale_,2*numberRows_,
     1194                           2*newNumberRows,-COIN_DBL_MAX,false);
     1195  savedColumnScale_ = resizeDouble(savedColumnScale_,2*numberColumns_,
     1196                           2*newNumberColumns,-COIN_DBL_MAX,false);
    10381197  if (objective_)
    10391198    objective_->resize(newNumberColumns);
     
    10591218    ray_ = NULL;
    10601219  }
    1061   delete [] rowScale_;
    1062   rowScale_ = NULL;
    1063   delete [] columnScale_;
    1064   columnScale_ = NULL;
     1220  setRowScale(NULL);
     1221  setColumnScale(NULL);
    10651222  if (status_) {
    10661223    if (newNumberColumns+newNumberRows) {
     
    11831340  delete [] ray_;
    11841341  ray_ = NULL;
    1185   delete [] rowScale_;
     1342  if (savedRowScale_!=rowScale_) {
     1343    delete [] rowScale_;
     1344    delete [] columnScale_;
     1345  }
    11861346  rowScale_ = NULL;
    1187   delete [] columnScale_;
    11881347  columnScale_ = NULL;
    11891348}
     
    12641423  delete [] ray_;
    12651424  ray_ = NULL;
    1266   delete [] rowScale_;
    1267   rowScale_ = NULL;
    1268   delete [] columnScale_;
    1269   columnScale_ = NULL;
     1425  setRowScale(NULL);
     1426  setColumnScale(NULL);
    12701427}
    12711428// Add one row
     
    13231480    if (!matrix_)
    13241481      createEmptyMatrix();
    1325     delete [] rowScale_;
    1326     rowScale_ = NULL;
    1327     delete [] columnScale_;
    1328     columnScale_ = NULL;
     1482    setRowScale(NULL);
     1483    setColumnScale(NULL);
    13291484#ifndef CLP_NO_STD
    13301485    if (lengthNames_) {
     
    14151570  if (rows)
    14161571    matrix_->appendRows(number,rows);
    1417   delete [] rowScale_;
    1418   rowScale_ = NULL;
    1419   delete [] columnScale_;
    1420   columnScale_ = NULL;
     1572  setRowScale(NULL);
     1573  setColumnScale(NULL);
    14211574  if (lengthNames_) {
    14221575    rowNames_.resize(numberRows_);
     
    18041957    if (!matrix_)
    18051958      createEmptyMatrix();
    1806     delete [] rowScale_;
    1807     rowScale_ = NULL;
    1808     delete [] columnScale_;
    1809     columnScale_ = NULL;
     1959    setRowScale(NULL);
     1960    setColumnScale(NULL);
    18101961#ifndef CLP_NO_STD
    18111962    if (lengthNames_) {
     
    19082059  if (columns)
    19092060    matrix_->appendCols(number,columns);
    1910   delete [] rowScale_;
    1911   rowScale_ = NULL;
    1912   delete [] columnScale_;
    1913   columnScale_ = NULL;
     2061  setRowScale(NULL);
     2062  setColumnScale(NULL);
    19142063  if (lengthNames_) {
    19152064    columnNames_.resize(numberColumns_);
     
    27912940                     int numberColumns, const int * whichColumn,
    27922941                     bool dropNames, bool dropIntegers)
     2942  :  specialOptions_(rhs->specialOptions_),
     2943  maximumColumns_(-1),
     2944  maximumRows_(-1),
     2945  savedRowScale_(NULL),
     2946  savedColumnScale_(NULL)
    27932947{
    27942948  defaultHandler_ = rhs->defaultHandler_;
     
    28012955  messages_ = rhs->messages_;
    28022956  coinMessages_ = rhs->coinMessages_;
     2957  maximumColumns_ = -1;
     2958  maximumRows_ = -1;
     2959  savedRowScale_ = NULL;
     2960  savedColumnScale_ = NULL;
    28032961  intParam_[ClpMaxNumIteration] = rhs->intParam_[ClpMaxNumIteration];
    28042962  intParam_[ClpMaxNumIterationHotStart] =
     
    33943552  } else if (!mode) {
    33953553    scalingFlag_=0;
    3396     delete [] rowScale_;
    3397     rowScale_ = NULL;
    3398     delete [] columnScale_;
    3399     columnScale_ = NULL;
     3554    setRowScale(NULL);
     3555    setColumnScale(NULL);
    34003556  }
    34013557}
     
    34763632 
    34773633  scalingFlag_=0;
    3478   delete [] rowScale_;
    3479   rowScale_ = NULL;
    3480   delete [] columnScale_;
    3481   columnScale_ = NULL;
     3634  setRowScale(NULL);
     3635  setColumnScale(NULL);
    34823636}
    34833637void
     
    35823736  }
    35833737  return coinModel;
     3738}
     3739// Start or reset using maximumRows_ and Columns_
     3740void
     3741ClpModel::startPermanentArrays()
     3742{
     3743  if ((specialOptions_&65536)!=0) {
     3744    if (numberRows_>maximumRows_||numberColumns_>maximumColumns_) {
     3745      if (numberRows_>maximumRows_)
     3746        maximumRows_ = numberRows_+10+numberRows_/100;
     3747      if (numberColumns_>maximumColumns_)
     3748        maximumColumns_ = numberColumns_+10+numberColumns_/100;
     3749      // need to make sure numberRows_ OK and size of matrices
     3750      resize(maximumRows_,maximumColumns_);
     3751    } else {
     3752      return;
     3753    }
     3754  } else {
     3755    specialOptions_ |= 65536;
     3756    maximumRows_ = numberRows_;
     3757    maximumColumns_ = numberColumns_;
     3758  }
     3759}
     3760// Stop using maximumRows_ and Columns_
     3761void
     3762ClpModel::stopPermanentArrays()
     3763{
     3764  specialOptions_ &= ~65536;
     3765  maximumRows_=-1;
     3766  maximumColumns_=-1;
     3767  if (rowScale_!=savedRowScale_) {
     3768    delete [] savedRowScale_;
     3769    delete [] savedColumnScale_;
     3770  }
     3771  savedRowScale_ = NULL;
     3772  savedColumnScale_ = NULL;
    35843773}
    35853774/* Find a network subset.
  • trunk/Clp/src/ClpModel.hpp

    r1141 r1147  
    258258  */
    259259  int cleanMatrix(double threshold=1.0e-20);
     260  /// Copy contents - resizing if necessary - otherwise re-use memory
     261  void copy(const ClpMatrixBase * from, ClpMatrixBase * & to);
    260262#ifndef CLP_NO_STD
    261263  /// Drops names - makes lengthnames 0 and names empty
     
    517519   inline const double * rowScale() const {return rowScale_;}
    518520   inline const double * columnScale() const {return columnScale_;}
    519    inline void setRowScale(double * scale) { delete [] (double *) rowScale_; rowScale_ = scale;}
    520    inline void setColumnScale(double * scale) { delete [] (double *) columnScale_; columnScale_ = scale;}
     521   inline double * mutableRowScale() const {return rowScale_;}
     522   inline double * mutableColumnScale() const {return columnScale_;}
     523   void setRowScale(double * scale) ;
     524   void setColumnScale(double * scale);
    521525  /// Scaling of objective
    522526  inline double objectiveScale() const
     
    827831      8192 - Do Primal when cleaning up primal
    828832      16384 - In fast dual (so we can switch off things)
    829       32678 - called from Osi
    830       65356 - keep arrays around as much as possible
     833      32768 - called from Osi
     834      65536 - keep arrays around as much as possible (also use maximumR/C)
    831835      131072 - scale factor arrays have inverse values at end
    832836      NOTE - many applications can call Clp but there may be some short cuts
     
    849853  //@{
    850854protected:
    851   /// Does most of deletion
    852   void gutsOfDelete();
     855  /// Does most of deletion (0 = all, 1 = most)
     856  void gutsOfDelete(int type);
    853857  /** Does most of copying
    854       If trueCopy false then just points to arrays */
    855   void gutsOfCopy(const ClpModel & rhs, bool trueCopy=true);
     858      If trueCopy 0 then just points to arrays
     859      If -1 leaves as much as possible */
     860  void gutsOfCopy(const ClpModel & rhs, int trueCopy=1);
    856861  /// gets lower and upper bounds on rows
    857862  void getRowBound(int iRow, double& lower, double& upper) const;
     
    864869  /// Does much of scaling
    865870  void gutsOfScaling();
    866    /// Objective value - always minimize
    867    inline double rawObjectiveValue() const {
    868       return objectiveValue_;
    869    }
     871  /// Objective value - always minimize
     872  inline double rawObjectiveValue() const {
     873    return objectiveValue_;
     874  }
     875  /// If we are using maximumRows_ and Columns_
     876  inline bool permanentArrays() const
     877  { return (specialOptions_&65536)!=0;}
     878  /// Start using maximumRows_ and Columns_
     879  void startPermanentArrays();
     880  /// Stop using maximumRows_ and Columns_
     881  void stopPermanentArrays();
    870882  /// Create row names as char **
    871883  const char * const * rowNamesAsChar() const;
     
    9961008  /// Coin messages
    9971009  CoinMessages coinMessages_;
     1010  /// Maximum number of columns in model
     1011  int maximumColumns_;
     1012  /// Maximum number of rows in model
     1013  int maximumRows_;
     1014  /// Base packed matrix
     1015  CoinPackedMatrix baseMatrix_;
     1016  /// Base row copy
     1017  CoinPackedMatrix baseRowCopy_;
     1018  /// Saved row scale factors for matrix
     1019  double * savedRowScale_;
     1020  /// Saved column scale factors
     1021  double * savedColumnScale_;
    9981022#ifndef CLP_NO_STD
    9991023  /// Array of string parameters
  • trunk/Clp/src/ClpPackedMatrix.cpp

    r1141 r1147  
    152152{
    153153  return new ClpPackedMatrix(*this);
     154}
     155// Copy contents - resizing if necessary - otherwise re-use memory
     156void
     157ClpPackedMatrix::copy(const ClpPackedMatrix * rhs)
     158{
     159  //*this = *rhs;
     160  assert(numberActiveColumns_ == rhs->numberActiveColumns_);
     161  assert (matrix_->isColOrdered()==rhs->matrix_->isColOrdered());
     162  matrix_->copyReuseArrays(*rhs->matrix_);
    154163}
    155164/* Subset clone (without gaps).  Duplicates are allowed
     
    20852094  int numberColumns = matrix_->getNumCols();
    20862095  // If empty - return as sanityCheck will trap
    2087   if (!numberRows||!numberColumns)
     2096  if (!numberRows||!numberColumns) {
     2097    model->setRowScale(NULL);
     2098    model->setColumnScale(NULL);
    20882099    return 1;
     2100  }
    20892101  ClpMatrixBase * rowCopyBase=model->rowCopy();
    20902102  if (!rowCopyBase) {
     
    21062118  const double * element = rowCopy->getElements();
    21072119  int scaleLength = ((model->specialOptions()&131072)==0) ? 1 : 2;
    2108   double * rowScale = new double [numberRows*scaleLength];
    2109   double * columnScale = new double [numberColumns*scaleLength];
     2120  double * rowScale;
     2121  double * columnScale;
     2122  //assert (!model->rowScale());
     2123  bool arraysExist;
     2124  if (!model->rowScale()) {
     2125    rowScale = new double [numberRows*scaleLength];
     2126    columnScale = new double [numberColumns*scaleLength];
     2127    arraysExist=false;
     2128  } else {
     2129    rowScale=model->mutableRowScale();
     2130    columnScale=model->mutableColumnScale();
     2131    arraysExist=true;
     2132  }
    21102133  // we are going to mark bits we are interested in
    21112134  char * usefulRow = new char [numberRows];
     
    21642187    model->messageHandler()->message(CLP_PACKEDSCALE_FORGET,*model->messagesPointer())
    21652188      <<CoinMessageEol;
    2166     delete [] rowScale;
     2189    if (!arraysExist) {
     2190      delete [] rowScale;
     2191      delete [] columnScale;
     2192    } else {
     2193      model->setRowScale(NULL);
     2194      model->setColumnScale(NULL);
     2195    }
    21672196    delete [] usefulRow;
    2168     delete [] columnScale;
    21692197    delete [] usefulColumn;
    21702198    if (!model->rowCopy())
     
    25172545        inverseScale[iColumn] = 1.0/columnScale[iColumn] ;
    25182546    }
    2519     model->setRowScale(rowScale);
    2520     model->setColumnScale(columnScale);
     2547    if (!arraysExist) {
     2548      model->setRowScale(rowScale);
     2549      model->setColumnScale(columnScale);
     2550    }
    25212551    if (model->rowCopy()) {
    25222552      // need to replace row by row
  • trunk/Clp/src/ClpPackedMatrix.hpp

    r1055 r1147  
    303303  /// Clone
    304304  virtual ClpMatrixBase * clone() const ;
     305  /// Copy contents - resizing if necessary - otherwise re-use memory
     306  virtual void copy(const ClpPackedMatrix * from);
    305307  /** Subset clone (without gaps).  Duplicates are allowed
    306308      and order is as given */
  • trunk/Clp/src/ClpSimplex.cpp

    r1141 r1147  
    119119  allowedInfeasibility_(10.0),
    120120  automaticScale_(0),
    121   progress_(NULL)
     121  baseModel_(NULL)
    122122{
    123123  int i;
     
    145145// Subproblem constructor
    146146ClpSimplex::ClpSimplex ( const ClpModel * rhs,
    147                     int numberRows, const int * whichRow,
    148                     int numberColumns, const int * whichColumn,
    149                     bool dropNames, bool dropIntegers,bool fixOthers)
     147                        int numberRows, const int * whichRow,
     148                        int numberColumns, const int * whichColumn,
     149                        bool dropNames, bool dropIntegers,bool fixOthers)
    150150  : ClpModel(rhs, numberRows, whichRow,
    151151             numberColumns,whichColumn,dropNames,dropIntegers),
    152   columnPrimalInfeasibility_(0.0),
    153   rowPrimalInfeasibility_(0.0),
    154   columnPrimalSequence_(-2),
    155   rowPrimalSequence_(-2),
    156   columnDualInfeasibility_(0.0),
    157   rowDualInfeasibility_(0.0),
    158   moreSpecialOptions_(2),
    159   baseIteration_(0),
    160   primalToleranceToGetOptimal_(-1.0),
    161   remainingDualInfeasibility_(0.0),
    162   largeValue_(1.0e15),
    163   largestPrimalError_(0.0),
    164   largestDualError_(0.0),
    165   alphaAccuracy_(-1.0),
    166   dualBound_(1.0e10),
    167   alpha_(0.0),
    168   theta_(0.0),
     152    columnPrimalInfeasibility_(0.0),
     153    rowPrimalInfeasibility_(0.0),
     154    columnPrimalSequence_(-2),
     155    rowPrimalSequence_(-2),
     156    columnDualInfeasibility_(0.0),
     157    rowDualInfeasibility_(0.0),
     158    moreSpecialOptions_(2),
     159    baseIteration_(0),
     160    primalToleranceToGetOptimal_(-1.0),
     161    remainingDualInfeasibility_(0.0),
     162    largeValue_(1.0e15),
     163    largestPrimalError_(0.0),
     164    largestDualError_(0.0),
     165    alphaAccuracy_(-1.0),
     166    dualBound_(1.0e10),
     167    alpha_(0.0),
     168    theta_(0.0),
    169169  lowerIn_(0.0),
    170170  valueIn_(0.0),
     
    211211  pivotVariable_(NULL),
    212212  factorization_(NULL),
    213   savedSolution_(NULL),
     213    savedSolution_(NULL),
    214214  numberTimesOptimal_(0),
    215215  disasterArea_(NULL),
     
    229229  incomingInfeasibility_(1.0),
    230230  allowedInfeasibility_(10.0),
    231   automaticScale_(0),
    232   progress_(NULL)
     231    automaticScale_(0),
     232  baseModel_(NULL)
    233233{
    234234  int i;
     
    15251525  int out = sequenceOut_;
    15261526  matrix_->correctSequence(this,in,out);
    1527   int cycle=progress_->cycle(in,out,
     1527  int cycle=progress_.cycle(in,out,
    15281528                            directionIn_,directionOut_);
    15291529  if (cycle>0&&objective_->type()<2) {
     
    15321532      printf("Cycle of %d\n",cycle);
    15331533    // reset
    1534     progress_->startCheck();
     1534    progress_.startCheck();
    15351535    double random = randomNumberGenerator_.randomDouble();
    15361536    int extra = (int) (9.999*random);
     
    15581558  // only time to re-factorize if one before real time
    15591559  // this is so user won't be surprised that maximumPivots has exact meaning
    1560   CoinBigIndex lengthL = factorization_->numberElementsL();
    1561   CoinBigIndex lengthU = factorization_->numberElementsU();
    1562   CoinBigIndex lengthR = factorization_->numberElementsR();
    15631560  int numberPivots=factorization_->pivots();
    15641561  int maximumPivots = factorization_->maximumPivots();
     
    15781575    }
    15791576    return 1;
    1580   } else if (numberPivots*3>maximumPivots*2&&
    1581              lengthR*3>(lengthL+lengthU)*2+1000&&!numberDense) {
     1577  } else if (factorization_->timeToRefactorize()) {
    15821578    //printf("ret after %d pivots\n",factorization_->pivots());
    15831579    return 1;
     
    16891685  allowedInfeasibility_(10.0),
    16901686  automaticScale_(0),
    1691   progress_(NULL)
     1687  baseModel_(NULL)
    16921688{
    16931689  int i;
     
    17931789  allowedInfeasibility_(10.0),
    17941790  automaticScale_(0),
    1795   progress_(NULL)
     1791  baseModel_(NULL)
    17961792{
    17971793  int i;
     
    18361832  int numberRows2 = numberRows_+numberExtraRows_;
    18371833  auxiliaryModel_ = NULL;
     1834  moreSpecialOptions_ = rhs.moreSpecialOptions_;
    18381835  if ((whatsChanged_&1)!=0) {
    1839     lower_ = ClpCopyOfArray(rhs.lower_,numberColumns_+numberRows2);
     1836    int numberTotal = numberColumns_+numberRows2;
     1837    if ((specialOptions_&65536)!=0&&maximumRows_>=0) {
     1838      assert (maximumRows_>=numberRows2);
     1839      assert (maximumColumns_>=numberColumns_);
     1840      numberTotal = 2*(maximumColumns_+ maximumRows_);
     1841    }
     1842    lower_ = ClpCopyOfArray(rhs.lower_,numberTotal);
    18401843    rowLowerWork_ = lower_+numberColumns_;
    18411844    columnLowerWork_ = lower_;
    1842     upper_ = ClpCopyOfArray(rhs.upper_,numberColumns_+numberRows2);
     1845    upper_ = ClpCopyOfArray(rhs.upper_,numberTotal);
    18431846    rowUpperWork_ = upper_+numberColumns_;
    18441847    columnUpperWork_ = upper_;
    1845     //cost_ = ClpCopyOfArray(rhs.cost_,2*(numberColumns_+numberRows_));
    1846     cost_ = ClpCopyOfArray(rhs.cost_,(numberColumns_+numberRows2));
     1848    cost_ = ClpCopyOfArray(rhs.cost_,numberTotal);
    18471849    objectiveWork_ = cost_;
    18481850    rowObjectiveWork_ = cost_+numberColumns_;
    1849     dj_ = ClpCopyOfArray(rhs.dj_,numberRows2+numberColumns_);
     1851    dj_ = ClpCopyOfArray(rhs.dj_,numberTotal);
    18501852    if (dj_) {
    18511853      reducedCostWork_ = dj_;
    18521854      rowReducedCost_ = dj_+numberColumns_;
    18531855    }
    1854     solution_ = ClpCopyOfArray(rhs.solution_,numberRows2+numberColumns_);
     1856    solution_ = ClpCopyOfArray(rhs.solution_,numberTotal);
    18551857    if (solution_) {
    18561858      columnActivityWork_ = solution_;
     
    19121914  rowPrimalSequence_ = rhs.rowPrimalSequence_;
    19131915  columnDualInfeasibility_ = rhs.columnDualInfeasibility_;
    1914   moreSpecialOptions_ = rhs.moreSpecialOptions_;
    19151916  rowDualInfeasibility_ = rhs.rowDualInfeasibility_;
    19161917  baseIteration_ = rhs.baseIteration_;
     
    19651966  allowedInfeasibility_ = rhs.allowedInfeasibility_;
    19661967  automaticScale_ = rhs.automaticScale_;
    1967   if (rhs.progress_)
    1968     progress_ = new ClpSimplexProgress(*rhs.progress_);
    1969   else
    1970     progress_=NULL;
     1968  if (rhs.baseModel_) {
     1969    baseModel_ = new ClpSimplex(*rhs.baseModel_);
     1970  } else {
     1971    baseModel_ = NULL;
     1972  }
     1973  progress_=rhs.progress_;
    19711974  for (int i=0;i<4;i++) {
    19721975    spareIntArray_[i]=rhs.spareIntArray_[i];
     
    19861989ClpSimplex::gutsOfDelete(int type)
    19871990{
    1988   delete [] lower_;
    1989   lower_=NULL;
    1990   rowLowerWork_=NULL;
    1991   columnLowerWork_=NULL;
    1992   delete [] upper_;
    1993   upper_=NULL;
    1994   rowUpperWork_=NULL;
    1995   columnUpperWork_=NULL;
    1996   delete [] cost_;
    1997   cost_=NULL;
    1998   objectiveWork_=NULL;
    1999   rowObjectiveWork_=NULL;
    2000   delete [] dj_;
    2001   dj_=NULL;
    2002   reducedCostWork_=NULL;
    2003   rowReducedCost_=NULL;
    2004   delete [] solution_;
    2005   solution_=NULL;
    2006   rowActivityWork_=NULL;
    2007   columnActivityWork_=NULL;
    2008   delete [] savedSolution_;
    2009   savedSolution_ = NULL;
     1991  if (!type||(specialOptions_&65536)==0) {
     1992    delete [] lower_;
     1993    lower_=NULL;
     1994    rowLowerWork_=NULL;
     1995    columnLowerWork_=NULL;
     1996    delete [] upper_;
     1997    upper_=NULL;
     1998    rowUpperWork_=NULL;
     1999    columnUpperWork_=NULL;
     2000    delete [] cost_;
     2001    cost_=NULL;
     2002    objectiveWork_=NULL;
     2003    rowObjectiveWork_=NULL;
     2004    delete [] dj_;
     2005    dj_=NULL;
     2006    reducedCostWork_=NULL;
     2007    rowReducedCost_=NULL;
     2008    delete [] solution_;
     2009    solution_=NULL;
     2010    rowActivityWork_=NULL;
     2011    columnActivityWork_=NULL;
     2012    delete [] savedSolution_;
     2013    savedSolution_ = NULL;
     2014  }
    20102015  if ((specialOptions_&2)==0) {
    20112016    delete nonLinearCost_;
     
    20362041    delete primalColumnPivot_;
    20372042    primalColumnPivot_ = NULL;
    2038     delete progress_;
    2039     progress_=NULL;
     2043    delete baseModel_;
     2044    baseModel_=NULL;
    20402045  } else {
    20412046    // delete any size information in methods
     
    22702275    firstFree_ = firstFreeDual;
    22712276  } else if (numberSuperBasicWithDj||
    2272              (progress_&&progress_->lastIterationNumber(0)<=0)) {
     2277             (progress_.lastIterationNumber(0)<=0)) {
    22732278    firstFree_=firstFreePrimal;
    22742279  }
     
    23872392    firstFree_ = firstFreeDual;
    23882393  } else if (numberSuperBasicWithDj||
    2389              (progress_&&progress_->lastIterationNumber(0)<=0)) {
     2394             (progress_.lastIterationNumber(0)<=0)) {
    23902395    firstFree_=firstFreePrimal;
    23912396  }
     
    24972502  }
    24982503  if (what==63) {
     2504    pivotRow_=-1;
    24992505    if (!status_)
    25002506      createStatus();
     
    25422548  int numberRows2 = numberRows_+numberExtraRows_;
    25432549  int numberTotal = numberRows2+numberColumns_;
     2550  if ((specialOptions_&65536)!=0&&!auxiliaryModel_) {
     2551    assert (!numberExtraRows_);
     2552    if (!cost_||numberRows2>maximumRows_||numberColumns_>maximumColumns_) {
     2553      newArrays=true;
     2554      keepPivots=false;
     2555      if (cost_) {
     2556        if (numberRows2>maximumRows_)
     2557          maximumRows_ = numberRows2+10+numberRows2/100;
     2558        if (numberColumns_>maximumColumns_)
     2559          maximumColumns_ = numberColumns_+10+numberColumns_/100;
     2560      } else {
     2561        maximumRows_ = numberRows2;
     2562        maximumColumns_ = numberColumns_;
     2563      }
     2564      int numberTotal2 = (maximumRows_+maximumColumns_)*2;
     2565      delete [] cost_;
     2566      cost_ = new double[numberTotal2];
     2567      delete [] lower_;
     2568      delete [] upper_;
     2569      lower_ = new double[numberTotal2];
     2570      upper_ = new double[numberTotal2];
     2571      delete [] dj_;
     2572      dj_ = new double[numberTotal2];
     2573      delete [] solution_;
     2574      solution_ = new double[numberTotal2];
     2575      // ***** should be non NULL but seems to be too much
     2576      //printf("resize %d savedRowScale %x\n",maximumRows_,savedRowScale_);
     2577      if (savedRowScale_) {
     2578        delete [] savedRowScale_;
     2579        savedRowScale_ = new double [2*maximumRows_];
     2580        delete [] savedColumnScale_;
     2581        savedColumnScale_ = new double [2*maximumColumns_];
     2582      }
     2583    }
     2584  }
    25442585  int i;
    25452586  bool doSanityCheck=true;
     
    25942635    }
    25952636    if (scalingFlag_>0&&!rowScale_) {
     2637      if ((specialOptions_&65536)!=0&&!auxiliaryModel_) {
     2638        assert (!rowScale_);
     2639        rowScale_ = savedRowScale_;
     2640        columnScale_ = savedColumnScale_;
     2641      }
    25962642      if (matrix_->scale(this))
    25972643        scalingFlag_=-scalingFlag_; // not scaled after all
     
    27242770  }
    27252771  if (what==63) {
    2726     if (newArrays) {
     2772    if (newArrays&&(specialOptions_&65536)==0) {
    27272773      delete [] cost_;
    2728       // extra copy with original costs
    2729       //cost_ = new double[2*numberTotal];
    27302774      cost_ = new double[numberTotal];
    27312775      delete [] lower_;
    27322776      delete [] upper_;
    2733       lower_ = new double[numberColumns_+numberRows2];
    2734       upper_ = new double[numberColumns_+numberRows2];
     2777      lower_ = new double[numberTotal];
     2778      upper_ = new double[numberTotal];
    27352779      delete [] dj_;
    2736       dj_ = new double[numberRows2+numberColumns_];
     2780      dj_ = new double[numberTotal];
    27372781      delete [] solution_;
    2738       solution_ = new double[numberRows2+numberColumns_];
     2782      solution_ = new double[numberTotal];
    27392783    } else if (auxiliaryModel_) {
    27402784      assert (!cost_);
     
    35353579  if ((what&5)!=0)
    35363580    matrix_->generalExpanded(this,9,what); // update costs and bounds if necessary
     3581  if (goodMatrix&&(specialOptions_&65536)!=0&&!auxiliaryModel_) {
     3582    int save = maximumColumns_+maximumRows_;
     3583    memcpy(cost_+save,cost_,numberTotal*sizeof(double));
     3584    memcpy(lower_+save,lower_,numberTotal*sizeof(double));
     3585    memcpy(upper_+save,upper_,numberTotal*sizeof(double));
     3586    memcpy(dj_+save,dj_,numberTotal*sizeof(double));
     3587    memcpy(solution_+save,solution_,numberTotal*sizeof(double));
     3588    if (rowScale_&&false) {
     3589      memcpy(savedRowScale_+maximumRows_,rowScale_,numberRows2*sizeof(double));
     3590      memcpy(savedColumnScale_+maximumColumns_,columnScale_,numberColumns_*sizeof(double));
     3591    }
     3592  }
    35373593  return goodMatrix;
     3594}
     3595// Does rows and columns
     3596void
     3597ClpSimplex::createRim1(bool initial)
     3598{
     3599  int i;
     3600  int numberRows2 = numberRows_+numberExtraRows_;
     3601  int numberTotal = numberRows2+numberColumns_;
     3602  //assert (!auxiliaryModel_);
     3603  if (!auxiliaryModel_) {
     3604    if ((specialOptions_&65536)!=0&&true) {
     3605      assert (!initial);
     3606      int save = maximumColumns_+maximumRows_;
     3607      memcpy(lower_,lower_+save,numberTotal*sizeof(double));
     3608      memcpy(upper_,upper_+save,numberTotal*sizeof(double));
     3609      return;
     3610    }
     3611    const double * rowScale = rowScale_;
     3612    const double * columnScale = columnScale_;
     3613    // clean up any mismatches on infinity
     3614    // and fix any variables with tiny gaps
     3615    double primalTolerance=dblParam_[ClpPrimalTolerance];
     3616    if(rowScale) {
     3617      // If scaled then do all columns later in one loop
     3618      if (!initial) {
     3619        if ((specialOptions_&131072)==0) {
     3620          for (i=0;i<numberColumns_;i++) {
     3621            double multiplier = rhsScale_/columnScale[i];
     3622            double lowerValue = columnLower_[i];
     3623            double upperValue = columnUpper_[i];
     3624            if (lowerValue>-1.0e20) {
     3625              columnLowerWork_[i]=lowerValue*multiplier;
     3626              if (upperValue>=1.0e20) {
     3627                columnUpperWork_[i]=COIN_DBL_MAX;
     3628              } else {
     3629                columnUpperWork_[i]=upperValue*multiplier;
     3630                if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
     3631                  if (columnLowerWork_[i]>=0.0) {
     3632                    columnUpperWork_[i] = columnLowerWork_[i];
     3633                  } else if (columnUpperWork_[i]<=0.0) {
     3634                    columnLowerWork_[i] = columnUpperWork_[i];
     3635                  } else {
     3636                    columnUpperWork_[i] = 0.0;
     3637                    columnLowerWork_[i] = 0.0;
     3638                  }
     3639                }
     3640              }
     3641            } else if (upperValue<1.0e20) {
     3642              columnLowerWork_[i]=-COIN_DBL_MAX;
     3643              columnUpperWork_[i]=upperValue*multiplier;
     3644            } else {
     3645              // free
     3646              columnLowerWork_[i]=-COIN_DBL_MAX;
     3647              columnUpperWork_[i]=COIN_DBL_MAX;
     3648            }
     3649          }
     3650        } else {
     3651          const double * inverseScale = columnScale_+numberColumns_;
     3652          for (i=0;i<numberColumns_;i++) {
     3653            double multiplier = rhsScale_*inverseScale[i];
     3654            double lowerValue = columnLower_[i];
     3655            double upperValue = columnUpper_[i];
     3656            if (lowerValue>-1.0e20) {
     3657              columnLowerWork_[i]=lowerValue*multiplier;
     3658              if (upperValue>=1.0e20) {
     3659                columnUpperWork_[i]=COIN_DBL_MAX;
     3660              } else {
     3661                columnUpperWork_[i]=upperValue*multiplier;
     3662                if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
     3663                  if (columnLowerWork_[i]>=0.0) {
     3664                    columnUpperWork_[i] = columnLowerWork_[i];
     3665                  } else if (columnUpperWork_[i]<=0.0) {
     3666                    columnLowerWork_[i] = columnUpperWork_[i];
     3667                  } else {
     3668                    columnUpperWork_[i] = 0.0;
     3669                    columnLowerWork_[i] = 0.0;
     3670                  }
     3671                }
     3672              }
     3673            } else if (upperValue<1.0e20) {
     3674              columnLowerWork_[i]=-COIN_DBL_MAX;
     3675              columnUpperWork_[i]=upperValue*multiplier;
     3676            } else {
     3677              // free
     3678              columnLowerWork_[i]=-COIN_DBL_MAX;
     3679              columnUpperWork_[i]=COIN_DBL_MAX;
     3680            }
     3681          }
     3682        }
     3683      }
     3684      for (i=0;i<numberRows_;i++) {
     3685        double multiplier = rhsScale_*rowScale[i];
     3686        double lowerValue = rowLower_[i];
     3687        double upperValue = rowUpper_[i];
     3688        if (lowerValue>-1.0e20) {
     3689          rowLowerWork_[i]=lowerValue*multiplier;
     3690          if (upperValue>=1.0e20) {
     3691            rowUpperWork_[i]=COIN_DBL_MAX;
     3692          } else {
     3693            rowUpperWork_[i]=upperValue*multiplier;
     3694            if (fabs(rowUpperWork_[i]-rowLowerWork_[i])<=primalTolerance) {
     3695              if (rowLowerWork_[i]>=0.0) {
     3696                rowUpperWork_[i] = rowLowerWork_[i];
     3697              } else if (rowUpperWork_[i]<=0.0) {
     3698                rowLowerWork_[i] = rowUpperWork_[i];
     3699              } else {
     3700                rowUpperWork_[i] = 0.0;
     3701                rowLowerWork_[i] = 0.0;
     3702              }
     3703            }
     3704          }
     3705        } else if (upperValue<1.0e20) {
     3706          rowLowerWork_[i]=-COIN_DBL_MAX;
     3707          rowUpperWork_[i]=upperValue*multiplier;
     3708        } else {
     3709          // free
     3710          rowLowerWork_[i]=-COIN_DBL_MAX;
     3711          rowUpperWork_[i]=COIN_DBL_MAX;
     3712        }
     3713      }
     3714    } else if (rhsScale_!=1.0) {
     3715      for (i=0;i<numberColumns_;i++) {
     3716        double lowerValue = columnLower_[i];
     3717        double upperValue = columnUpper_[i];
     3718        if (lowerValue>-1.0e20) {
     3719          columnLowerWork_[i]=lowerValue*rhsScale_;
     3720          if (upperValue>=1.0e20) {
     3721            columnUpperWork_[i]=COIN_DBL_MAX;
     3722          } else {
     3723            columnUpperWork_[i]=upperValue*rhsScale_;
     3724            if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
     3725              if (columnLowerWork_[i]>=0.0) {
     3726                columnUpperWork_[i] = columnLowerWork_[i];
     3727              } else if (columnUpperWork_[i]<=0.0) {
     3728                columnLowerWork_[i] = columnUpperWork_[i];
     3729              } else {
     3730                columnUpperWork_[i] = 0.0;
     3731                columnLowerWork_[i] = 0.0;
     3732              }
     3733            }
     3734          }
     3735        } else if (upperValue<1.0e20) {
     3736          columnLowerWork_[i]=-COIN_DBL_MAX;
     3737          columnUpperWork_[i]=upperValue*rhsScale_;
     3738        } else {
     3739          // free
     3740          columnLowerWork_[i]=-COIN_DBL_MAX;
     3741          columnUpperWork_[i]=COIN_DBL_MAX;
     3742        }
     3743      }
     3744      for (i=0;i<numberRows_;i++) {
     3745        double lowerValue = rowLower_[i];
     3746        double upperValue = rowUpper_[i];
     3747        if (lowerValue>-1.0e20) {
     3748          rowLowerWork_[i]=lowerValue*rhsScale_;
     3749          if (upperValue>=1.0e20) {
     3750            rowUpperWork_[i]=COIN_DBL_MAX;
     3751          } else {
     3752            rowUpperWork_[i]=upperValue*rhsScale_;
     3753            if (fabs(rowUpperWork_[i]-rowLowerWork_[i])<=primalTolerance) {
     3754              if (rowLowerWork_[i]>=0.0) {
     3755                rowUpperWork_[i] = rowLowerWork_[i];
     3756              } else if (rowUpperWork_[i]<=0.0) {
     3757                rowLowerWork_[i] = rowUpperWork_[i];
     3758              } else {
     3759                rowUpperWork_[i] = 0.0;
     3760                rowLowerWork_[i] = 0.0;
     3761              }
     3762            }
     3763          }
     3764        } else if (upperValue<1.0e20) {
     3765          rowLowerWork_[i]=-COIN_DBL_MAX;
     3766          rowUpperWork_[i]=upperValue*rhsScale_;
     3767        } else {
     3768          // free
     3769          rowLowerWork_[i]=-COIN_DBL_MAX;
     3770          rowUpperWork_[i]=COIN_DBL_MAX;
     3771        }
     3772      }
     3773    } else {
     3774      for (i=0;i<numberColumns_;i++) {
     3775        double lowerValue = columnLower_[i];
     3776        double upperValue = columnUpper_[i];
     3777        if (lowerValue>-1.0e20) {
     3778          columnLowerWork_[i]=lowerValue;
     3779          if (upperValue>=1.0e20) {
     3780            columnUpperWork_[i]=COIN_DBL_MAX;
     3781          } else {
     3782            columnUpperWork_[i]=upperValue;
     3783            if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
     3784              if (columnLowerWork_[i]>=0.0) {
     3785                columnUpperWork_[i] = columnLowerWork_[i];
     3786              } else if (columnUpperWork_[i]<=0.0) {
     3787                columnLowerWork_[i] = columnUpperWork_[i];
     3788              } else {
     3789                columnUpperWork_[i] = 0.0;
     3790                columnLowerWork_[i] = 0.0;
     3791              }
     3792            }
     3793          }
     3794        } else if (upperValue<1.0e20) {
     3795          columnLowerWork_[i]=-COIN_DBL_MAX;
     3796          columnUpperWork_[i]=upperValue;
     3797        } else {
     3798          // free
     3799          columnLowerWork_[i]=-COIN_DBL_MAX;
     3800          columnUpperWork_[i]=COIN_DBL_MAX;
     3801        }
     3802      }
     3803      for (i=0;i<numberRows_;i++) {
     3804        double lowerValue = rowLower_[i];
     3805        double upperValue = rowUpper_[i];
     3806        if (lowerValue>-1.0e20) {
     3807          rowLowerWork_[i]=lowerValue;
     3808          if (upperValue>=1.0e20) {
     3809            rowUpperWork_[i]=COIN_DBL_MAX;
     3810          } else {
     3811            rowUpperWork_[i]=upperValue;
     3812            if (fabs(rowUpperWork_[i]-rowLowerWork_[i])<=primalTolerance) {
     3813              if (rowLowerWork_[i]>=0.0) {
     3814                rowUpperWork_[i] = rowLowerWork_[i];
     3815              } else if (rowUpperWork_[i]<=0.0) {
     3816                rowLowerWork_[i] = rowUpperWork_[i];
     3817              } else {
     3818                rowUpperWork_[i] = 0.0;
     3819                rowLowerWork_[i] = 0.0;
     3820              }
     3821            }
     3822          }
     3823        } else if (upperValue<1.0e20) {
     3824          rowLowerWork_[i]=-COIN_DBL_MAX;
     3825          rowUpperWork_[i]=upperValue;
     3826        } else {
     3827          // free
     3828          rowLowerWork_[i]=-COIN_DBL_MAX;
     3829          rowUpperWork_[i]=COIN_DBL_MAX;
     3830        }
     3831      }
     3832    }
     3833  } else {
     3834    // auxiliary model
     3835    if (!initial) {
     3836      // just copy
     3837      memcpy(lower_,auxiliaryModel_->lower_+numberTotal,numberTotal*sizeof(double));
     3838      memcpy(upper_,auxiliaryModel_->upper_+numberTotal,numberTotal*sizeof(double));
     3839    } else {
     3840      // clean up any mismatches on infinity
     3841      // and fix any variables with tiny gaps
     3842      const double * rowScale = auxiliaryModel_->rowScale_ ;
     3843      const double * columnScale = auxiliaryModel_->columnScale_ ;
     3844      assert (rhsScale_==1.0);
     3845      assert (objectiveScale_==1.0);
     3846      if ((auxiliaryModel_->specialOptions_&2)==0) {
     3847        // need to do column bounds
     3848        if (columnScale) {
     3849          const double * inverseScale = columnScale+numberColumns_;
     3850          // scaled
     3851          for (i=0;i<numberColumns_;i++) {
     3852            double value;
     3853            value = columnLower_[i];
     3854            if (value>-1.0e20)
     3855              value *= inverseScale[i];
     3856            lower_[i] = value;
     3857            value = columnUpper_[i];
     3858            if (value<1.0e20)
     3859              value *= inverseScale[i];
     3860            upper_[i] = value;
     3861          }
     3862        } else {
     3863          for (i=0;i<numberColumns_;i++) {
     3864            lower_[i] = columnLower_[i];
     3865            upper_[i] = columnUpper_[i];
     3866          }
     3867        }
     3868        // save
     3869        memcpy(auxiliaryModel_->lower_+numberTotal,lower_,numberColumns_*sizeof(double));
     3870        memcpy(auxiliaryModel_->upper_+numberTotal,upper_,numberColumns_*sizeof(double));
     3871      } else {
     3872        // just copy
     3873        memcpy(lower_,auxiliaryModel_->lower_+numberTotal,numberColumns_*sizeof(double));
     3874        memcpy(upper_,auxiliaryModel_->upper_+numberTotal,numberColumns_*sizeof(double));
     3875      }
     3876      if ((auxiliaryModel_->specialOptions_&1)==0) {
     3877        // need to do row bounds
     3878        if (rowScale) {
     3879          // scaled
     3880          for (i=0;i<numberRows_;i++) {
     3881            double value;
     3882            value = rowLower_[i];
     3883            if (value>-1.0e20)
     3884              value *= rowScale[i];
     3885            lower_[i+numberColumns_] = value;
     3886            value = rowUpper_[i];
     3887            if (value<1.0e20)
     3888              value *= rowScale[i];
     3889            upper_[i+numberColumns_] = value;
     3890          }
     3891        } else {
     3892          for (i=0;i<numberRows_;i++) {
     3893            lower_[i+numberColumns_] = rowLower_[i];
     3894            upper_[i+numberColumns_] = rowUpper_[i];
     3895          }
     3896        }
     3897        // save
     3898        memcpy(auxiliaryModel_->lower_+numberTotal+numberColumns_,
     3899               lower_+numberColumns_,numberRows_*sizeof(double));
     3900        memcpy(auxiliaryModel_->upper_+numberTotal+numberColumns_,
     3901               upper_+numberColumns_,numberRows_*sizeof(double));
     3902      } else {
     3903        // just copy
     3904        memcpy(lower_+numberColumns_,
     3905               auxiliaryModel_->lower_+numberTotal+numberColumns_,numberRows_*sizeof(double));
     3906        memcpy(upper_+numberColumns_,
     3907               auxiliaryModel_->upper_+numberTotal+numberColumns_,numberRows_*sizeof(double));
     3908      }
     3909    }
     3910    if (initial) {
     3911      // do basis
     3912      assert ((auxiliaryModel_->specialOptions_&8)!=0);
     3913      // clean solution trusting basis
     3914      for ( i=0;i<numberTotal;i++) {
     3915        Status status =getStatus(i);
     3916        double value=solution_[i]; // may as well leave if basic (values pass)
     3917        if (status!=basic) {
     3918          if (upper_[i]==lower_[i]) {
     3919            setStatus(i,isFixed);
     3920            value=lower_[i];
     3921          } else if (status==atLowerBound) {
     3922            if (lower_[i]>-1.0e20) {
     3923              value=lower_[i];
     3924            } else {
     3925              if (upper_[i]<1.0e20) {
     3926                value=upper_[i];
     3927                setStatus(i,atUpperBound);
     3928              } else {
     3929                setStatus(i,isFree);
     3930              }
     3931            }
     3932          } else if (status==atUpperBound) {
     3933            if (upper_[i]<1.0e20) {
     3934              value=upper_[i];
     3935            } else {
     3936              if (lower_[i]>-1.0e20) {
     3937                value=lower_[i];
     3938                setStatus(i,atLowerBound);
     3939              } else {
     3940                setStatus(i,isFree);
     3941              }
     3942            }
     3943          }
     3944        }
     3945        solution_[i]=value;
     3946      }
     3947    }
     3948  }
     3949  if ((specialOptions_&65536)!=0&&false) {
     3950    assert (!initial);
     3951    int save = maximumColumns_+maximumRows_;
     3952    for (int i=0;i<numberTotal;i++) {
     3953      assert (fabs(lower_[i]-lower_[i+save])<1.0e-5);
     3954      assert (fabs(upper_[i]-upper_[i+save])<1.0e-5);
     3955    }
     3956  }
     3957}
     3958// Does objective
     3959void
     3960ClpSimplex::createRim4(bool initial)
     3961{
     3962  int i;
     3963  int numberRows2 = numberRows_+numberExtraRows_;
     3964  int numberTotal = numberRows2+numberColumns_;
     3965  //assert (!auxiliaryModel_);
     3966  if (!auxiliaryModel_||(initial&&(auxiliaryModel_->specialOptions_&4)==0)) {
     3967    if ((specialOptions_&65536)!=0&&true) {
     3968      assert (!initial);
     3969      int save = maximumColumns_+maximumRows_;
     3970      memcpy(cost_,cost_+save,numberTotal*sizeof(double));
     3971      return;
     3972    }
     3973    double direction = optimizationDirection_*objectiveScale_;
     3974    const double * obj = objective();
     3975    const double * rowScale = auxiliaryModel_ ? auxiliaryModel_->rowScale_ : rowScale_;
     3976    const double * columnScale = auxiliaryModel_ ? auxiliaryModel_->columnScale_ : columnScale_;
     3977    // and also scale by scale factors
     3978    if (rowScale) {
     3979      if (rowObjective_) {
     3980        for (i=0;i<numberRows_;i++)
     3981          rowObjectiveWork_[i] = rowObjective_[i]*direction/rowScale[i];
     3982      } else {
     3983        memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
     3984      }
     3985      // If scaled then do all columns later in one loop
     3986      if (!initial||auxiliaryModel_) {
     3987        for (i=0;i<numberColumns_;i++) {
     3988          CoinAssert(fabs(obj[i])<1.0e25);
     3989          objectiveWork_[i] = obj[i]*direction*columnScale[i];
     3990        }
     3991      }
     3992    } else {
     3993      if (rowObjective_) {
     3994        for (i=0;i<numberRows_;i++)
     3995          rowObjectiveWork_[i] = rowObjective_[i]*direction;
     3996      } else {
     3997        memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
     3998      }
     3999      for (i=0;i<numberColumns_;i++) {
     4000        CoinAssert(fabs(obj[i])<1.0e25);
     4001        objectiveWork_[i] = obj[i]*direction;
     4002      }
     4003    }
     4004    if (auxiliaryModel_) {
     4005      // save costs
     4006      memcpy(auxiliaryModel_->cost_+numberTotal,cost_,numberTotal*sizeof(double));
     4007    }
     4008  } else {
     4009    // just copy
     4010    memcpy(cost_,auxiliaryModel_->cost_+numberTotal,numberTotal*sizeof(double));
     4011  }
     4012}
     4013// Does rows and columns and objective
     4014void
     4015ClpSimplex::createRim5(bool initial)
     4016{
     4017  createRim4(initial);
     4018  createRim1(initial);
    35384019}
    35394020void
     
    52595740  primalColumnPivot_ = otherModel.primalColumnPivot_->clone(true);
    52605741  perturbation_ = otherModel.perturbation_;
    5261   specialOptions_ = otherModel.specialOptions_;
    52625742  moreSpecialOptions_ = otherModel.moreSpecialOptions_;
    52635743  automaticScale_ = otherModel.automaticScale_;
     
    55185998    // save event handler in case already set
    55195999    ClpEventHandler * handler = eventHandler_->clone();
    5520     ClpModel::gutsOfDelete();
     6000    ClpModel::gutsOfDelete(0);
    55216001    eventHandler_ = handler;
    55226002    gutsOfDelete(0);
     
    71847664    if (!useFactorization||factorization_->numberRows()!=numberRows_) {
    71857665      useFactorization = false;
    7186       factorization_->increasingRows(2);
    7187       // row activities have negative sign
    7188       factorization_->slackValue(-1.0);
    7189       factorization_->zeroTolerance(1.0e-13);
     7666      factorization_->setDefaultValues();
    71907667      // Switch off dense (unless special option set)
    71917668      if ((specialOptions_&8)==0)
     
    73117788  saved.objectiveScale_ = objectiveScale_;
    73127789  // Progress indicator
    7313   delete progress_;
    7314   progress_ = new ClpSimplexProgress (this);
     7790  progress_.fillFromModel (this);
    73157791  return saved;
    73167792}
     
    73277803  objectiveScale_ = saved.objectiveScale_;
    73287804  acceptablePivot_ = saved.acceptablePivot_;
    7329   delete progress_;
    7330   progress_=NULL;
    73317805}
    73327806// To flag a variable (not inline to allow for column generation)
     
    74887962  CoinAssert(incomingInfeasibility_>=0.0);
    74897963  CoinAssert(allowedInfeasibility_>=incomingInfeasibility_);
    7490 }
    7491 //#############################################################################
    7492 
    7493 ClpSimplexProgress::ClpSimplexProgress ()
    7494 {
    7495   int i;
    7496   for (i=0;i<CLP_PROGRESS;i++) {
    7497     objective_[i] = COIN_DBL_MAX;
    7498     infeasibility_[i] = -1.0; // set to an impossible value
    7499     realInfeasibility_[i] = COIN_DBL_MAX;
    7500     numberInfeasibilities_[i]=-1;
    7501     iterationNumber_[i]=-1;
    7502   }
    7503 #ifdef CLP_PROGRESS_WEIGHT
    7504   for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
    7505     objectiveWeight_[i] = COIN_DBL_MAX;
    7506     infeasibilityWeight_[i] = -1.0; // set to an impossible value
    7507     realInfeasibilityWeight_[i] = COIN_DBL_MAX;
    7508     numberInfeasibilitiesWeight_[i]=-1;
    7509     iterationNumberWeight_[i]=-1;
    7510   }
    7511   drop_ =0.0;
    7512   best_ =0.0;
    7513 #endif
    7514   initialWeight_=0.0;
    7515   for (i=0;i<CLP_CYCLE;i++) {
    7516     //obj_[i]=COIN_DBL_MAX;
    7517     in_[i]=-1;
    7518     out_[i]=-1;
    7519     way_[i]=0;
    7520   }
    7521   numberTimes_ = 0;
    7522   numberBadTimes_ = 0;
    7523   model_ = NULL;
    7524   oddState_=0;
    7525 }
    7526 
    7527 
    7528 //-----------------------------------------------------------------------------
    7529 
    7530 ClpSimplexProgress::~ClpSimplexProgress ()
    7531 {
    7532 }
    7533 // Copy constructor.
    7534 ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs)
    7535 {
    7536   int i;
    7537   for (i=0;i<CLP_PROGRESS;i++) {
    7538     objective_[i] = rhs.objective_[i];
    7539     infeasibility_[i] = rhs.infeasibility_[i];
    7540     realInfeasibility_[i] = rhs.realInfeasibility_[i];
    7541     numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i];
    7542     iterationNumber_[i]=rhs.iterationNumber_[i];
    7543   }
    7544 #ifdef CLP_PROGRESS_WEIGHT
    7545   for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
    7546     objectiveWeight_[i] = rhs.objectiveWeight_[i];
    7547     infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
    7548     realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
    7549     numberInfeasibilitiesWeight_[i]=rhs.numberInfeasibilitiesWeight_[i];
    7550     iterationNumberWeight_[i]=rhs.iterationNumberWeight_[i];
    7551   }
    7552   drop_ = rhs.drop_;
    7553   best_ = rhs.best_;
    7554 #endif
    7555   initialWeight_ = rhs.initialWeight_;
    7556   for (i=0;i<CLP_CYCLE;i++) {
    7557     //obj_[i]=rhs.obj_[i];
    7558     in_[i]=rhs.in_[i];
    7559     out_[i]=rhs.out_[i];
    7560     way_[i]=rhs.way_[i];
    7561   }
    7562   numberTimes_ = rhs.numberTimes_;
    7563   numberBadTimes_ = rhs.numberBadTimes_;
    7564   model_ = rhs.model_;
    7565   oddState_=rhs.oddState_;
    7566 }
    7567 // Copy constructor.from model
    7568 ClpSimplexProgress::ClpSimplexProgress(ClpSimplex * model)
    7569 {
    7570   model_ = model;
    7571   reset();
    7572   initialWeight_=0.0;
    7573 }
    7574 // Assignment operator. This copies the data
    7575 ClpSimplexProgress &
    7576 ClpSimplexProgress::operator=(const ClpSimplexProgress & rhs)
    7577 {
    7578   if (this != &rhs) {
    7579     int i;
    7580     for (i=0;i<CLP_PROGRESS;i++) {
    7581       objective_[i] = rhs.objective_[i];
    7582       infeasibility_[i] = rhs.infeasibility_[i];
    7583       realInfeasibility_[i] = rhs.realInfeasibility_[i];
    7584       numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i];
    7585       iterationNumber_[i]=rhs.iterationNumber_[i];
    7586     }
    7587 #ifdef CLP_PROGRESS_WEIGHT
    7588     for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
    7589       objectiveWeight_[i] = rhs.objectiveWeight_[i];
    7590       infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
    7591       realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
    7592       numberInfeasibilitiesWeight_[i]=rhs.numberInfeasibilitiesWeight_[i];
    7593       iterationNumberWeight_[i]=rhs.iterationNumberWeight_[i];
    7594     }
    7595     drop_ = rhs.drop_;
    7596     best_ = rhs.best_;
    7597 #endif
    7598     initialWeight_ = rhs.initialWeight_;
    7599     for (i=0;i<CLP_CYCLE;i++) {
    7600       //obj_[i]=rhs.obj_[i];
    7601       in_[i]=rhs.in_[i];
    7602       out_[i]=rhs.out_[i];
    7603       way_[i]=rhs.way_[i];
    7604     }
    7605     numberTimes_ = rhs.numberTimes_;
    7606     numberBadTimes_ = rhs.numberBadTimes_;
    7607     model_ = rhs.model_;
    7608     oddState_=rhs.oddState_;
    7609   }
    7610   return *this;
    7611 }
    7612 // Seems to be something odd about exact comparison of doubles on linux
    7613 static bool equalDouble(double value1, double value2)
    7614 {
    7615 
    7616   union { double d; int i[2]; } v1,v2;
    7617   v1.d = value1;
    7618   v2.d = value2;
    7619   if (sizeof(int)*2==sizeof(double))
    7620     return (v1.i[0]==v2.i[0]&&v1.i[1]==v2.i[1]);
    7621   else
    7622     return (v1.i[0]==v2.i[0]);
    7623 }
    7624 int
    7625 ClpSimplexProgress::looping()
    7626 {
    7627   if (!model_)
    7628     return -1;
    7629   double objective = model_->rawObjectiveValue();
    7630   double infeasibility;
    7631   double realInfeasibility=0.0;
    7632   int numberInfeasibilities;
    7633   int iterationNumber = model_->numberIterations();
    7634   if (model_->algorithm()<0) {
    7635     // dual
    7636     infeasibility = model_->sumPrimalInfeasibilities();
    7637     numberInfeasibilities = model_->numberPrimalInfeasibilities();
    7638   } else {
    7639     //primal
    7640     infeasibility = model_->sumDualInfeasibilities();
    7641     realInfeasibility = model_->nonLinearCost()->sumInfeasibilities();
    7642     numberInfeasibilities = model_->numberDualInfeasibilities();
    7643   }
    7644   int i;
    7645   int numberMatched=0;
    7646   int matched=0;
    7647   int nsame=0;
    7648   for (i=0;i<CLP_PROGRESS;i++) {
    7649     bool matchedOnObjective = equalDouble(objective,objective_[i]);
    7650     bool matchedOnInfeasibility = equalDouble(infeasibility,infeasibility_[i]);
    7651     bool matchedOnInfeasibilities =
    7652       (numberInfeasibilities==numberInfeasibilities_[i]);
    7653    
    7654     if (matchedOnObjective&&matchedOnInfeasibility&&matchedOnInfeasibilities) {
    7655       matched |= (1<<i);
    7656       // Check not same iteration
    7657       if (iterationNumber!=iterationNumber_[i]) {
    7658         numberMatched++;
    7659         // here mainly to get over compiler bug?
    7660         if (model_->messageHandler()->logLevel()>10)
    7661           printf("%d %d %d %d %d loop check\n",i,numberMatched,
    7662                  matchedOnObjective, matchedOnInfeasibility,
    7663                  matchedOnInfeasibilities);
    7664       } else {
    7665         // stuck but code should notice
    7666         nsame++;
    7667       }
    7668     }
    7669     if (i) {
    7670       objective_[i-1] = objective_[i];
    7671       infeasibility_[i-1] = infeasibility_[i];
    7672       realInfeasibility_[i-1] = realInfeasibility_[i];
    7673       numberInfeasibilities_[i-1]=numberInfeasibilities_[i];
    7674       iterationNumber_[i-1]=iterationNumber_[i];
    7675     }
    7676   }
    7677   objective_[CLP_PROGRESS-1] = objective;
    7678   infeasibility_[CLP_PROGRESS-1] = infeasibility;
    7679   realInfeasibility_[CLP_PROGRESS-1] = realInfeasibility;
    7680   numberInfeasibilities_[CLP_PROGRESS-1]=numberInfeasibilities;
    7681   iterationNumber_[CLP_PROGRESS-1]=iterationNumber;
    7682   if (nsame==CLP_PROGRESS)
    7683     numberMatched=CLP_PROGRESS; // really stuck
    7684   if (model_->progressFlag())
    7685     numberMatched=0;
    7686   numberTimes_++;
    7687   if (numberTimes_<10)
    7688     numberMatched=0;
    7689   // skip if just last time as may be checking something
    7690   if (matched==(1<<(CLP_PROGRESS-1)))
    7691     numberMatched=0;
    7692   if (numberMatched) {
    7693     model_->messageHandler()->message(CLP_POSSIBLELOOP,model_->messages())
    7694       <<numberMatched
    7695       <<matched
    7696       <<numberTimes_
    7697       <<CoinMessageEol;
    7698     numberBadTimes_++;
    7699     if (numberBadTimes_<10) {
    7700       // make factorize every iteration
    7701       model_->forceFactorization(1);
    7702       if (numberBadTimes_<2) {
    7703         startCheck(); // clear other loop check
    7704         if (model_->algorithm()<0) {
    7705           // dual - change tolerance
    7706           model_->setCurrentDualTolerance(model_->currentDualTolerance()*1.05);
    7707           // if infeasible increase dual bound
    7708           if (model_->dualBound()<1.0e17) {
    7709             model_->setDualBound(model_->dualBound()*1.1);
    7710           }
    7711         } else {
    7712           // primal - change tolerance 
    7713           if (numberBadTimes_>3)
    7714             model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance()*1.05);
    7715           // if infeasible increase infeasibility cost
    7716           if (model_->nonLinearCost()->numberInfeasibilities()&&
    7717               model_->infeasibilityCost()<1.0e17) {
    7718             model_->setInfeasibilityCost(model_->infeasibilityCost()*1.1);
    7719           }
    7720         }
    7721       } else {
    7722         // flag
    7723         int iSequence;
    7724         if (model_->algorithm()<0) {
    7725           // dual
    7726           if (model_->dualBound()>1.0e14)
    7727             model_->setDualBound(1.0e14);
    7728           iSequence=in_[CLP_CYCLE-1];
    7729         } else {
    7730           // primal
    7731           if (model_->infeasibilityCost()>1.0e14)
    7732             model_->setInfeasibilityCost(1.0e14);
    7733           iSequence=out_[CLP_CYCLE-1];
    7734         }
    7735         if (iSequence>=0) {
    7736           char x = model_->isColumn(iSequence) ? 'C' :'R';
    7737           if (model_->messageHandler()->logLevel()>=63)
    7738             model_->messageHandler()->message(CLP_SIMPLEX_FLAG,model_->messages())
    7739               <<x<<model_->sequenceWithin(iSequence)
    7740               <<CoinMessageEol;
    7741           // if Gub then needs to be sequenceIn_
    7742           int save=model_->sequenceIn();
    7743           model_->setSequenceIn(iSequence);
    7744           model_->setFlagged(iSequence);
    7745           model_->setSequenceIn(save);
    7746           //printf("flagging %d from loop\n",iSequence);
    7747           startCheck();
    7748         } else {
    7749           // Give up
    7750           if (model_->messageHandler()->logLevel()>=63)
    7751             printf("***** All flagged?\n");
    7752           return 4;
    7753         }
    7754         // reset
    7755         numberBadTimes_=2;
    7756       }
    7757       return -2;
    7758     } else {
    7759       // look at solution and maybe declare victory
    7760       if (infeasibility<1.0e-4) {
    7761         return 0;
    7762       } else {
    7763         model_->messageHandler()->message(CLP_LOOP,model_->messages())
    7764           <<CoinMessageEol;
    7765 #ifndef NDEBUG
    7766         printf("debug loop ClpSimplex A\n");
    7767         abort();
    7768 #endif
    7769         return 3;
    7770       }
    7771     }
    7772   }
    7773   return -1;
    7774 }
    7775 // Resets as much as possible
    7776 void
    7777 ClpSimplexProgress::reset()
    7778 {
    7779   int i;
    7780   for (i=0;i<CLP_PROGRESS;i++) {
    7781     if (model_->algorithm()>=0)
    7782       objective_[i] = COIN_DBL_MAX;
    7783     else
    7784       objective_[i] = -COIN_DBL_MAX;
    7785     infeasibility_[i] = -1.0; // set to an impossible value
    7786     realInfeasibility_[i] = COIN_DBL_MAX;
    7787     numberInfeasibilities_[i]=-1;
    7788     iterationNumber_[i]=-1;
    7789   }
    7790 #ifdef CLP_PROGRESS_WEIGHT
    7791   for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
    7792     objectiveWeight_[i] = COIN_DBL_MAX;
    7793     infeasibilityWeight_[i] = -1.0; // set to an impossible value
    7794     realInfeasibilityWeight_[i] = COIN_DBL_MAX;
    7795     numberInfeasibilitiesWeight_[i]=-1;
    7796     iterationNumberWeight_[i]=-1;
    7797   }
    7798   drop_ =0.0;
    7799   best_ =0.0;
    7800 #endif
    7801   for (i=0;i<CLP_CYCLE;i++) {
    7802     //obj_[i]=COIN_DBL_MAX;
    7803     in_[i]=-1;
    7804     out_[i]=-1;
    7805     way_[i]=0;
    7806   }
    7807   numberTimes_ = 0;
    7808   numberBadTimes_ = 0;
    7809   oddState_=0;
    7810 }
    7811 // Returns previous objective (if -1) - current if (0)
    7812 double
    7813 ClpSimplexProgress::lastObjective(int back) const
    7814 {
    7815   return objective_[CLP_PROGRESS-1-back];
    7816 }
    7817 // Returns previous infeasibility (if -1) - current if (0)
    7818 double
    7819 ClpSimplexProgress::lastInfeasibility(int back) const
    7820 {
    7821   return realInfeasibility_[CLP_PROGRESS-1-back];
    7822 }
    7823 // Sets real primal infeasibility
    7824 void
    7825 ClpSimplexProgress::setInfeasibility(double value)
    7826 {
    7827   for (int i=1;i<CLP_PROGRESS;i++)
    7828     realInfeasibility_[i-1] = realInfeasibility_[i];
    7829   realInfeasibility_[CLP_PROGRESS-1]=value;
    7830 }
    7831 // Modify objective e.g. if dual infeasible in dual
    7832 void
    7833 ClpSimplexProgress::modifyObjective(double value)
    7834 {
    7835   objective_[CLP_PROGRESS-1]=value;
    7836 }
    7837 // Returns previous iteration number (if -1) - current if (0)
    7838 int
    7839 ClpSimplexProgress::lastIterationNumber(int back) const
    7840 {
    7841   return iterationNumber_[CLP_PROGRESS-1-back];
    7842 }
    7843 // clears iteration numbers (to switch off panic)
    7844 void
    7845 ClpSimplexProgress::clearIterationNumbers()
    7846 {
    7847   for (int i=0;i<CLP_PROGRESS;i++)
    7848     iterationNumber_[i]=-1;
    7849 }
    7850 // Start check at beginning of whileIterating
    7851 void
    7852 ClpSimplexProgress::startCheck()
    7853 {
    7854   int i;
    7855   for (i=0;i<CLP_CYCLE;i++) {
    7856     //obj_[i]=COIN_DBL_MAX;
    7857     in_[i]=-1;
    7858     out_[i]=-1;
    7859     way_[i]=0;
    7860   }
    7861 }
    7862 // Returns cycle length in whileIterating
    7863 int
    7864 ClpSimplexProgress::cycle(int in, int out,int wayIn,int wayOut)
    7865 {
    7866   int i;
    7867 #if 0
    7868   if (model_->numberIterations()>206571) {
    7869     printf("in %d out %d\n",in,out);
    7870     for (i=0;i<CLP_CYCLE;i++)
    7871       printf("cy %d in %d out %d\n",i,in_[i],out_[i]);
    7872   }
    7873 #endif
    7874   int matched=0;
    7875   // first see if in matches any out
    7876   for (i=1;i<CLP_CYCLE;i++) {
    7877     if (in==out_[i]) {
    7878       // even if flip then suspicious
    7879       matched=-1;
    7880       break;
    7881     }
    7882   }
    7883 #if 0
    7884   if (!matched||in_[0]<0) {
    7885     // can't be cycle
    7886     for (i=0;i<CLP_CYCLE-1;i++) {
    7887       //obj_[i]=obj_[i+1];
    7888       in_[i]=in_[i+1];
    7889       out_[i]=out_[i+1];
    7890       way_[i]=way_[i+1];
    7891     }
    7892   } else {
    7893     // possible cycle
    7894     matched=0;
    7895     for (i=0;i<CLP_CYCLE-1;i++) {
    7896       int k;
    7897       char wayThis = way_[i];
    7898       int inThis = in_[i];
    7899       int outThis = out_[i];
    7900       //double objThis = obj_[i];
    7901       for(k=i+1;k<CLP_CYCLE;k++) {
    7902         if (inThis==in_[k]&&outThis==out_[k]&&wayThis==way_[k]) {
    7903           int distance = k-i;
    7904           if (k+distance<CLP_CYCLE) {
    7905             // See if repeats
    7906             int j=k+distance;
    7907             if (inThis==in_[j]&&outThis==out_[j]&&wayThis==way_[j]) {
    7908               matched=distance;
    7909               break;
    7910             }
    7911           } else {
    7912             matched=distance;
    7913             break;
    7914           }
    7915         }
    7916       }
    7917       //obj_[i]=obj_[i+1];
    7918       in_[i]=in_[i+1];
    7919       out_[i]=out_[i+1];
    7920       way_[i]=way_[i+1];
    7921     }
    7922   }
    7923 #else
    7924   if (matched&&in_[0]>=0) {
    7925     // possible cycle - only check [0] against all
    7926     matched=0;
    7927     int nMatched=0;
    7928     char way0 = way_[0];
    7929     int in0 = in_[0];
    7930     int out0 = out_[0];
    7931     //double obj0 = obj_[i];
    7932     for(int k=1;k<CLP_CYCLE-4;k++) {
    7933       if (in0==in_[k]&&out0==out_[k]&&way0==way_[k]) {
    7934         nMatched++;
    7935         // See if repeats
    7936         int end = CLP_CYCLE-k;
    7937         int j;
    7938         for ( j=1;j<end;j++) {
    7939           if (in_[j+k]!=in_[j]||out_[j+k]!=out_[j]||way_[j+k]!=way_[j])
    7940             break;
    7941         }
    7942         if (j==end) {
    7943           matched=k;
    7944           break;
    7945         }
    7946       }
    7947     }
    7948     // If three times then that is too much even if not regular
    7949     if (matched<=0&&nMatched>1)
    7950       matched=100;
    7951   }
    7952   for (i=0;i<CLP_CYCLE-1;i++) {
    7953     //obj_[i]=obj_[i+1];
    7954     in_[i]=in_[i+1];
    7955     out_[i]=out_[i+1];
    7956     way_[i]=way_[i+1];
    7957   }
    7958 #endif
    7959   char way = 1-wayIn+4*(1-wayOut);
    7960   //obj_[i]=model_->objectiveValue();
    7961   in_[CLP_CYCLE-1]=in;
    7962   out_[CLP_CYCLE-1]=out;
    7963   way_[CLP_CYCLE-1]=way;
    7964   return matched;
    79657964}
    79667965// Allow initial dense factorization
     
    87998798      // work arrays exist - update as well
    88008799      whatsChanged_ &= ~128;
     8800      double value;
    88018801      if (columnLower_[elementIndex]==-COIN_DBL_MAX) {
    8802         columnLowerWork_[elementIndex] = -COIN_DBL_MAX;
     8802        value = -COIN_DBL_MAX;
    88038803      } else if (!columnScale_) {
    8804         columnLowerWork_[elementIndex] = elementValue * rhsScale_;
     8804        value = elementValue * rhsScale_;
    88058805      } else {
    88068806        assert (!auxiliaryModel_);
    8807         columnLowerWork_[elementIndex] = elementValue * rhsScale_
     8807        value = elementValue * rhsScale_
    88088808          / columnScale_[elementIndex];
    88098809      }
     8810      columnLowerWork_[elementIndex] = value;
     8811      if (maximumRows_>=0)
     8812        columnLowerWork_[elementIndex+maximumRows_+maximumColumns_] = value;
    88108813    }
    88118814  }
     
    88308833      // work arrays exist - update as well
    88318834      whatsChanged_ &= ~256;
     8835      double value;
    88328836      if (columnUpper_[elementIndex]==COIN_DBL_MAX) {
    8833         columnUpperWork_[elementIndex] = COIN_DBL_MAX;
     8837        value = COIN_DBL_MAX;
    88348838      } else if (!columnScale_) {
    8835         columnUpperWork_[elementIndex] = elementValue * rhsScale_;
     8839        value = elementValue * rhsScale_;
    88368840      } else {
    88378841        assert (!auxiliaryModel_);
    8838         columnUpperWork_[elementIndex] = elementValue * rhsScale_
     8842        value = elementValue * rhsScale_
    88398843          / columnScale_[elementIndex];
    88408844      }
     8845      columnUpperWork_[elementIndex] = value;
     8846      if (maximumRows_>=0)
     8847        columnUpperWork_[elementIndex+maximumRows_+maximumColumns_] = value;
    88418848    }
    88428849  }
     
    93559362{
    93569363  if (value) {
    9357     specialOptions_|=65536;
     9364    //specialOptions_|=65536;
     9365    startPermanentArrays();
    93589366  } else {
    93599367    specialOptions_&=~65536;
     
    93819389    memcpy(rowActivity_,rhs.rowActivity_,numberRows_*sizeof(double));
    93829390    memcpy(dual_,rhs.dual_,numberRows_*sizeof(double));
     9391  }
     9392}
     9393// Save a copy of model with certain state - normally without cuts
     9394void
     9395ClpSimplex::makeBaseModel()
     9396{
     9397  delete baseModel_;
     9398  baseModel_ = new ClpSimplex(*this);
     9399}
     9400// Switch off base model
     9401void
     9402ClpSimplex::deleteBaseModel()
     9403{
     9404  delete baseModel_;
     9405  baseModel_=NULL;
     9406}
     9407// Reset to base model
     9408void
     9409ClpSimplex::setToBaseModel(ClpSimplex * model)
     9410{
     9411  if (!model)
     9412    model=baseModel_;
     9413  assert (model);
     9414  int multiplier = ((model->specialOptions_&65536)!=0) ? 2 : 1;
     9415  assert (multiplier==2);
     9416  if (multiplier==2) {
     9417    assert (model->maximumRows_>=0);
     9418    if (maximumRows_<0) {
     9419      specialOptions_ |= 65536;
     9420      maximumRows_ = model->maximumRows_;
     9421      maximumColumns_ = model->maximumColumns_;
     9422    }
     9423  }
     9424  // temporary - later use maximumRows_ for rowUpper_ etc
     9425  assert (numberRows_>=model->numberRows_);
     9426  abort();
     9427}
     9428// Start or reset using maximumRows_ and Columns_
     9429bool
     9430ClpSimplex::startPermanentArrays()
     9431{
     9432  int maximumRows = maximumRows_;
     9433  int maximumColumns = maximumColumns_;
     9434  ClpModel::startPermanentArrays();
     9435  if (maximumRows!=maximumRows_||
     9436      maximumColumns!=maximumColumns_) {
     9437    int numberTotal2 = (maximumRows_+maximumColumns_)*2;
     9438    delete [] cost_;
     9439    cost_ = new double[numberTotal2];
     9440    delete [] lower_;
     9441    delete [] upper_;
     9442    lower_ = new double[numberTotal2];
     9443    upper_ = new double[numberTotal2];
     9444    delete [] dj_;
     9445    dj_ = new double[numberTotal2];
     9446    delete [] solution_;
     9447    solution_ = new double[numberTotal2];
     9448    //delete [] savedRowScale_;
     9449    //savedRowScale_ = new double [2*maximumRows_];
     9450    //delete [] savedColumnScale_;
     9451    //savedColumnScale_ = new double [2*maximumColumns_];
     9452    return true;
     9453  } else {
     9454    return false;
    93839455  }
    93849456}
  • trunk/Clp/src/ClpSimplex.hpp

    r1137 r1147  
    2121class CoinIndexedVector;
    2222class ClpNonLinearCost;
    23 class ClpSimplexProgress;
    2423class CoinModel;
    2524class OsiClpSolverInterface;
     
    4342    There is an algorithm data member.  + for primal variations
    4443    and - for dual variations
    45 
    46     This file also includes (at end) a very simple class ClpSimplexProgress
    47     which is where anti-looping stuff should migrate to
    4844
    4945*/
     
    132128  inline bool usingAuxiliaryModel() const
    133129  { return auxiliaryModel_!=NULL;}
     130  /// Save a copy of model with certain state - normally without cuts
     131  void makeBaseModel();
     132  /// Switch off base model
     133  void deleteBaseModel();
     134  /// See if we have base model
     135  inline ClpSimplex *  baseModel() const
     136  { return baseModel_;}
     137  /** Reset to base model (just size and arrays needed)
     138      If model NULL use internal copy
     139  */
     140  void setToBaseModel(ClpSimplex * model=NULL);
    134141  /// Assignment operator. This copies the data
    135142    ClpSimplex & operator=(const ClpSimplex & rhs);
     
    735742  void gutsOfCopy(const ClpSimplex & rhs);
    736743  /** puts in format I like (rowLower,rowUpper) also see StandardMatrix
    737       1 bit does rows, 2 bit does column bounds, 4 bit does objective(s).
     744      1 bit does rows (now and columns), (2 bit does column bounds), 4 bit does objective(s).
    738745      8 bit does solution scaling in
    739746      16 bit does rowArray and columnArray indexed vectors
     
    746753  */
    747754  bool createRim(int what,bool makeRowCopy=false,int startFinishOptions=0);
     755  /// Does rows and columns
     756  void createRim1(bool initial);
     757  /// Does objective
     758  void createRim4(bool initial);
     759  /// Does rows and columns and objective
     760  void createRim5(bool initial);
    748761  /** releases above arrays and does solution scaling out.  May also
    749762      get rid of factorization data -
     
    788801    st_byte |= status;
    789802  }
     803  /// Start or reset using maximumRows_ and Columns_ - true if change
     804  bool startPermanentArrays();
    790805  /** Normally the first factorization does sparse coding because
    791806      the factorization could be singular.  This allows initial dense
     
    865880      1 bit - if presolve says infeasible in ClpSolve return
    866881      2 bit - if presolved problem infeasible return
     882      4 bit - keep arrays like upper_ around
    867883  */
    868884  inline int moreSpecialOptions() const
     
    871887      1 bit - if presolve says infeasible in ClpSolve return
    872888      2 bit - if presolved problem infeasible return
     889      4 bit - keep arrays like upper_ around
    873890  */
    874891  inline void setMoreSpecialOptions(int value)
     
    12751292  /// Automatic scaling of objective and rhs and bounds
    12761293  int automaticScale_;
     1294  /// A copy of model with certain state - normally without cuts
     1295  ClpSimplex * baseModel_;
    12771296  /// For dealing with all issues of cycling etc
    1278   ClpSimplexProgress * progress_;
     1297  ClpSimplexProgress progress_;
    12791298public:
    12801299  /// Spare int array for passing information [0]!=0 switches on
     
    12991318ClpSimplexUnitTest(const std::string & mpsDir);
    13001319
    1301 
    1302 /// For saving extra information to see if looping.
    1303 class ClpSimplexProgress {
    1304 
    1305 public:
    1306 
    1307 
    1308   /**@name Constructors and destructor and copy */
    1309   //@{
    1310   /// Default constructor
    1311     ClpSimplexProgress (  );
    1312 
    1313   /// Constructor from model
    1314     ClpSimplexProgress ( ClpSimplex * model );
    1315 
    1316   /// Copy constructor.
    1317   ClpSimplexProgress(const ClpSimplexProgress &);
    1318 
    1319   /// Assignment operator. This copies the data
    1320     ClpSimplexProgress & operator=(const ClpSimplexProgress & rhs);
    1321   /// Destructor
    1322    ~ClpSimplexProgress (  );
    1323   /// Resets as much as possible
    1324    void reset();
    1325   //@}
    1326 
    1327   /**@name Check progress */
    1328   //@{
    1329   /** Returns -1 if okay, -n+1 (n number of times bad) if bad but action taken,
    1330       >=0 if give up and use as problem status
    1331   */
    1332     int looping (  );
    1333   /// Start check at beginning of whileIterating
    1334   void startCheck();
    1335   /// Returns cycle length in whileIterating
    1336   int cycle(int in, int out,int wayIn,int wayOut);
    1337 
    1338   /// Returns previous objective (if -1) - current if (0)
    1339   double lastObjective(int back=1) const;
    1340   /// Set real primal infeasibility and move back
    1341   void setInfeasibility(double value);
    1342   /// Returns real primal infeasibility (if -1) - current if (0)
    1343   double lastInfeasibility(int back=1) const;
    1344   /// Modify objective e.g. if dual infeasible in dual
    1345   void modifyObjective(double value);
    1346   /// Returns previous iteration number (if -1) - current if (0)
    1347   int lastIterationNumber(int back=1) const;
    1348   /// clears all iteration numbers (to switch off panic)
    1349   void clearIterationNumbers();
    1350   /// Odd state
    1351   inline void newOddState()
    1352   { oddState_= - oddState_-1;}
    1353   inline void endOddState()
    1354   { oddState_=abs(oddState_);}
    1355   inline void clearOddState()
    1356   { oddState_=0;}
    1357   inline int oddState() const
    1358   { return oddState_;}
    1359   /// number of bad times
    1360   inline int badTimes() const
    1361   { return numberBadTimes_;}
    1362   inline void clearBadTimes()
    1363   { numberBadTimes_=0;}
    1364 
    1365   //@}
    1366   /**@name Data  */
    1367 #define CLP_PROGRESS 5
    1368   //#define CLP_PROGRESS_WEIGHT 10
    1369   //@{
    1370   /// Objective values
    1371   double objective_[CLP_PROGRESS];
    1372   /// Sum of infeasibilities for algorithm
    1373   double infeasibility_[CLP_PROGRESS];
    1374   /// Sum of real primal infeasibilities for primal
    1375   double realInfeasibility_[CLP_PROGRESS];
    1376 #ifdef CLP_PROGRESS_WEIGHT
    1377   /// Objective values for weights
    1378   double objectiveWeight_[CLP_PROGRESS_WEIGHT];
    1379   /// Sum of infeasibilities for algorithm for weights
    1380   double infeasibilityWeight_[CLP_PROGRESS_WEIGHT];
    1381   /// Sum of real primal infeasibilities for primal for weights
    1382   double realInfeasibilityWeight_[CLP_PROGRESS_WEIGHT];
    1383   /// Drop  for weights
    1384   double drop_;
    1385   /// Best? for weights
    1386   double best_;
    1387 #endif
    1388   /// Initial weight for weights
    1389   double initialWeight_;
    1390 #define CLP_CYCLE 12
    1391   /// For cycle checking
    1392   //double obj_[CLP_CYCLE];
    1393   int in_[CLP_CYCLE];
    1394   int out_[CLP_CYCLE];
    1395   char way_[CLP_CYCLE];
    1396   /// Pointer back to model so we can get information
    1397   ClpSimplex * model_;
    1398   /// Number of infeasibilities
    1399   int numberInfeasibilities_[CLP_PROGRESS];
    1400   /// Iteration number at which occurred
    1401   int iterationNumber_[CLP_PROGRESS];
    1402 #ifdef CLP_PROGRESS_WEIGHT
    1403   /// Number of infeasibilities for weights
    1404   int numberInfeasibilitiesWeight_[CLP_PROGRESS_WEIGHT];
    1405   /// Iteration number at which occurred for weights
    1406   int iterationNumberWeight_[CLP_PROGRESS_WEIGHT];
    1407 #endif
    1408   /// Number of times checked (so won't stop too early)
    1409   int numberTimes_;
    1410   /// Number of times it looked like loop
    1411   int numberBadTimes_;
    1412   /// If things are in an odd state
    1413   int oddState_;
    1414   //@}
    1415 };
    14161320// For Devex stuff
    14171321#define DEVEX_TRY_NORM 1.0e-4
  • trunk/Clp/src/ClpSimplexDual.cpp

    r1141 r1147  
    351351  int factorType=0;
    352352  // Start check for cycles
    353   progress_->startCheck();
     353  progress_.startCheck();
    354354  // Say change made on first iteration
    355355  changeMade_=1;
     
    429429    // If values pass then do easy ones on first time
    430430    if (ifValuesPass&&
    431         progress_->lastIterationNumber(0)<0&&saveDuals) {
     431        progress_.lastIterationNumber(0)<0&&saveDuals) {
    432432      doEasyOnesInValuesPass(saveDuals);
    433433    }
     
    625625    int factorType=0;
    626626    // Start check for cycles
    627     progress_->startCheck();
     627    progress_.startCheck();
    628628    // Say change made on first iteration
    629629    changeMade_=1;
     
    677677      // If values pass then do easy ones on first time
    678678      if (ifValuesPass&&
    679           progress_->lastIterationNumber(0)<0&&saveDuals) {
     679          progress_.lastIterationNumber(0)<0&&saveDuals) {
    680680        doEasyOnesInValuesPass(saveDuals);
    681681      }
     
    12041204#endif
    12051205              setFlagged(sequenceOut_);
    1206               progress_->clearBadTimes();
     1206              progress_.clearBadTimes();
    12071207              lastBadIteration_ = numberIterations_; // say be more cautious
    12081208              rowArray_[0]->clear();
     
    13291329#endif
    13301330            setFlagged(sequenceOut_);
    1331             progress_->clearBadTimes();
     1331            progress_.clearBadTimes();
    13321332            lastBadIteration_ = numberIterations_; // say be more cautious
    13331333            rowArray_[0]->clear();
     
    23262326    changeCost=0.0;
    23272327    // put back original bounds and then check
    2328     createRim(1);
     2328    createRim1(false);
    23292329    int iSequence;
    23302330    // bounds will get bigger - just look at ones at bounds
     
    34123412#endif
    34133413          setFlagged(sequenceOut_);
    3414           progress_->clearBadTimes();
     3414          progress_.clearBadTimes();
    34153415         
    34163416          // Go to safe
     
    34893489#endif
    34903490    setFlagged(sequenceOut_);
    3491     progress_->clearBadTimes();
     3491    progress_.clearBadTimes();
    34923492   
    34933493    // Go to safer
     
    35333533  }
    35343534  // Double check infeasibility if no action
    3535   if (progress_->lastIterationNumber(0)==numberIterations_) {
     3535  if (progress_.lastIterationNumber(0)==numberIterations_) {
    35363536    if (dualRowPivot_->looksOptimal()) {
    35373537      numberPrimalInfeasibilities_ = 0;
     
    35413541  } else {
    35423542    double thisObj = objectiveValue_;
    3543     double lastObj = progress_->lastObjective(0);
     3543    double lastObj = progress_.lastObjective(0);
    35443544    if(!ifValuesPass&&firstFree_<0) {
    35453545      if (lastObj>thisObj+1.0e-3*CoinMax(fabs(thisObj),fabs(lastObj))+1.0) {
     
    36063606  if (numberIterations_>CoinMax(1000,numberRows_>>4)&&(specialOptions_&64)!=0) {
    36073607    if (sumPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e5) {
    3608       int backIteration = progress_->lastIterationNumber(CLP_PROGRESS-1);
     3608      int backIteration = progress_.lastIterationNumber(CLP_PROGRESS-1);
    36093609      if (backIteration>0&&numberIterations_-backIteration<9*CLP_PROGRESS) {
    36103610        if (factorization_->pivotTolerance()<0.9) {
     
    36123612          factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance()*1.05+0.02,0.91));
    36133613          //printf("tol now %g\n",factorization_->pivotTolerance());
    3614           progress_->clearIterationNumbers();
     3614          progress_.clearIterationNumbers();
    36153615        }
    36163616      }
     
    36203620  int loop;
    36213621  if (!givenDuals&&type!=2)
    3622     loop = progress_->looping();
     3622    loop = progress_.looping();
    36233623  else
    36243624    loop=-1;
     
    36983698    //if (dualFeasible()||problemStatus_==-4||(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
    36993699    if (dualFeasible()||problemStatus_==-4) {
    3700       progress_->modifyObjective(objectiveValue_
     3700      progress_.modifyObjective(objectiveValue_
    37013701                               -sumDualInfeasibilities_*dualBound_);
    37023702      if (primalFeasible()&&!givenDuals) {
     
    37243724            changeBounds(true,NULL,changeCost);
    37253725            //computeObjectiveValue();
    3726             createRim(4);
     3726            createRim4(false);
    37273727            // make sure duals are current
    37283728            computeDuals(givenDuals);
     
    38073807              perturbation_=102; // stop any perturbations
    38083808              cleanDuals=1;
    3809               createRim(4);
     3809              createRim4(false);
    38103810              problemStatus_=-1;
    38113811            }
     
    38433843          // make sure fake bounds are back
    38443844          changeBounds(true,NULL,changeCost);
    3845           createRim(4);
     3845          createRim4(false);
    38463846        }
    38473847        if (numberChangedBounds<=0||dualBound_>1.0e20||
     
    38523852            //cleanDuals=1;
    38533853            //numberChangedBounds=1;
    3854             //createRim(4);
     3854            //createRim4(false);
    38553855          }
    38563856        } else {
     
    38893889        memcpy(xsolution,solution_,(numberRows_+numberColumns_)*sizeof(double));
    38903890#endif
    3891         createRim(4);
     3891        createRim4(false);
    38923892        // make sure duals are current
    38933893        computeDuals(givenDuals);
     
    41854185  }
    41864186#if 1
    4187   double thisObj = progress_->lastObjective(0);
    4188   double lastObj = progress_->lastObjective(1);
     4187  double thisObj = progress_.lastObjective(0);
     4188  double lastObj = progress_.lastObjective(1);
    41894189  if (lastObj>thisObj+1.0e-4*CoinMax(fabs(thisObj),fabs(lastObj))+1.0e-4
    41904190      &&givenDuals==NULL&&firstFree_<0) {
     
    47314731  if (!useFactorization||factorization_->numberRows()!=numberRows_) {
    47324732    useFactorization = false;
    4733     factorization_->increasingRows(2);
    4734     // row activities have negative sign
    4735     factorization_->slackValue(-1.0);
    4736     factorization_->zeroTolerance(1.0e-13);
     4733    factorization_->setDefaultValues();
    47374734
    47384735    int factorizationStatus = internalFactorize(0);
     
    51035100  }
    51045101  if (!useFactorization) {
    5105     factorization_->increasingRows(2);
    5106     // row activities have negative sign
    5107     factorization_->slackValue(-1.0);
    5108     factorization_->zeroTolerance(1.0e-13);
     5102    factorization_->setDefaultValues();
    51095103
    51105104    int factorizationStatus = internalFactorize(0);
     
    57975791ClpSimplexDual::resetFakeBounds()
    57985792{
     5793  // put back original bounds and then check
     5794  createRim1(false);
    57995795  double dummyChangeCost=0.0;
    58005796  changeBounds(true,rowArray_[2],dummyChangeCost);
  • trunk/Clp/src/ClpSimplexNonlinear.cpp

    r1141 r1147  
    7373    int factorType=0;
    7474    // Start check for cycles
    75     progress_->startCheck();
     75    progress_.startCheck();
    7676    /*
    7777      Status of problem:
     
    119119        }
    120120      }
    121       statusOfProblemInPrimal(lastCleaned,factorType,progress_,true,
     121      statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,
    122122                              bestObjectiveWhenFlagged);
    123123     
     
    25152515        <<CoinMessageEol;
    25162516      setFlagged(sequenceIn_);
    2517       progress_->clearBadTimes();
     2517      progress_.clearBadTimes();
    25182518      lastBadIteration_ = numberIterations_; // say be more cautious
    25192519      clearAll();
  • trunk/Clp/src/ClpSimplexOther.cpp

    r1128 r1147  
    19911991  int factorType=0;
    19921992  // Start check for cycles
    1993   progress_->startCheck();
     1993  progress_.startCheck();
    19941994  // Say change made on first iteration
    19951995  changeMade_=1;
     
    21072107  int loop;
    21082108  if (type!=2)
    2109     loop = progress_->looping();
     2109    loop = progress_.looping();
    21102110  else
    21112111    loop=-1;
     
    21562156  }
    21572157  if (dualFeasible()||problemStatus_==-4) {
    2158     progress_->modifyObjective(objectiveValue_
     2158    progress_.modifyObjective(objectiveValue_
    21592159                               -sumDualInfeasibilities_*dualBound_);
    21602160  }
     
    23272327                <<CoinMessageEol;
    23282328              setFlagged(sequenceOut_);
    2329               progress_->clearBadTimes();
     2329              progress_.clearBadTimes();
    23302330              lastBadIteration_ = numberIterations_; // say be more cautious
    23312331              rowArray_[0]->clear();
     
    24182418              <<CoinMessageEol;
    24192419            setFlagged(sequenceOut_);
    2420             progress_->clearBadTimes();
     2420            progress_.clearBadTimes();
    24212421            lastBadIteration_ = numberIterations_; // say be more cautious
    24222422            rowArray_[0]->clear();
     
    30043004  //double * rowLower = smallModel.rowLower();
    30053005  //double * rowUpper = smallModel.rowUpper();
    3006   const double * columnLower = smallModel.columnLower();
     3006  //const double * columnLower = smallModel.columnLower();
    30073007  //const double * columnUpper = smallModel.columnUpper();
    30083008  const CoinPackedMatrix * matrix = smallModel.matrix();
     
    30123012  const int * columnLength = matrix->getVectorLengths();
    30133013  const double * objective = smallModel.objective();
    3014   double objectiveOffset=0.0;
     3014  //double objectiveOffset=0.0;
    30153015  // would use for fixed?
    30163016  CoinZeroN(rhsOffset,numberRows_);
  • trunk/Clp/src/ClpSimplexPrimal.cpp

    r1141 r1147  
    257257    double lastObjectiveValue=COIN_DBL_MAX;
    258258    // Start check for cycles
    259     if (!progress_)
    260       progress_ = new ClpSimplexProgress(this);
    261     progress_->startCheck();
     259    progress_.fillFromModel(this);
     260    progress_.startCheck();
    262261    /*
    263262      Status of problem:
     
    311310         
    312311      // may factorize, checks if problem finished
    313       statusOfProblemInPrimal(lastCleaned,factorType,progress_,true,ifValuesPass,saveModel);
     312      statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,ifValuesPass,saveModel);
    314313      if (initialStatus==10) {
    315314        // cleanup phase
     
    320319              perturb(1);
    321320              matrix_->rhsOffset(this,true,false);
    322               statusOfProblemInPrimal(lastCleaned,factorType,progress_,true,ifValuesPass,saveModel);
     321              statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,ifValuesPass,saveModel);
    323322            }
    324323          }
     
    340339        //lastSprintIteration=numberIterations_;
    341340        printf("End small model after\n");
    342         statusOfProblemInPrimal(lastCleaned,factorType,progress_,true,ifValuesPass,saveModel);
     341        statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,ifValuesPass,saveModel);
    343342      }
    344343      int numberSprintIterations=0;
     
    350349        saveModel=NULL;
    351350        // Skip factorization
    352         //statusOfProblemInPrimal(lastCleaned,factorType,progress_,false,saveModel);
    353         statusOfProblemInPrimal(lastCleaned,factorType,progress_,true,ifValuesPass,saveModel);
     351        //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
     352        statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,ifValuesPass,saveModel);
    354353      } else if (problemStatus_<0&&!saveModel&&numberSprintColumns&&firstFree_<0) {
    355354        int numberSort=0;
     
    415414          delete [] weight;
    416415          // Skip factorization
    417           //statusOfProblemInPrimal(lastCleaned,factorType,progress_,false,saveModel);
    418           //statusOfProblemInPrimal(lastCleaned,factorType,progress_,true,saveModel);
     416          //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
     417          //statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,saveModel);
    419418          stopSprint = numberIterations_+numberSprintIterations;
    420419          printf("Sprint with %d columns for %d iterations\n",
     
    466465  // if infeasible get real values
    467466  //printf("XXXXY final cost %g\n",infeasibilityCost_);
    468   progress_->initialWeight_=0.0;
     467  progress_.initialWeight_=0.0;
    469468  if (problemStatus_==1&&secondaryStatus_!=6) {
    470469    infeasibilityCost_=0.0;
     
    26052604              <<CoinMessageEol;
    26062605            setFlagged(sequenceIn_);
    2607             progress_->clearBadTimes();
     2606            progress_.clearBadTimes();
    26082607            lastBadIteration_ = numberIterations_; // say be more cautious
    26092608            clearAll();
     
    26522651            <<CoinMessageEol;
    26532652          setFlagged(sequenceIn_);
    2654           progress_->clearBadTimes();
     2653          progress_.clearBadTimes();
    26552654          lastBadIteration_ = numberIterations_; // say be more cautious
    26562655          clearAll();
     
    27622761              <<CoinMessageEol;
    27632762            setFlagged(sequenceIn_);
    2764             progress_->clearBadTimes();
     2763            progress_.clearBadTimes();
    27652764          }
    27662765          lastBadIteration_ = numberIterations_; // say be more cautious
  • trunk/Clp/src/ClpSolve.cpp

    r1141 r1147  
    25842584  fprintf(fp,"3                    options,extraInfo,independentOptions);\n");
    25852585}
     2586//#############################################################################
     2587#include "ClpNonLinearCost.hpp"
     2588
     2589ClpSimplexProgress::ClpSimplexProgress ()
     2590{
     2591  int i;
     2592  for (i=0;i<CLP_PROGRESS;i++) {
     2593    objective_[i] = COIN_DBL_MAX;
     2594    infeasibility_[i] = -1.0; // set to an impossible value
     2595    realInfeasibility_[i] = COIN_DBL_MAX;
     2596    numberInfeasibilities_[i]=-1;
     2597    iterationNumber_[i]=-1;
     2598  }
     2599#ifdef CLP_PROGRESS_WEIGHT
     2600  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
     2601    objectiveWeight_[i] = COIN_DBL_MAX;
     2602    infeasibilityWeight_[i] = -1.0; // set to an impossible value
     2603    realInfeasibilityWeight_[i] = COIN_DBL_MAX;
     2604    numberInfeasibilitiesWeight_[i]=-1;
     2605    iterationNumberWeight_[i]=-1;
     2606  }
     2607  drop_ =0.0;
     2608  best_ =0.0;
     2609#endif
     2610  initialWeight_=0.0;
     2611  for (i=0;i<CLP_CYCLE;i++) {
     2612    //obj_[i]=COIN_DBL_MAX;
     2613    in_[i]=-1;
     2614    out_[i]=-1;
     2615    way_[i]=0;
     2616  }
     2617  numberTimes_ = 0;
     2618  numberBadTimes_ = 0;
     2619  model_ = NULL;
     2620  oddState_=0;
     2621}
     2622
     2623
     2624//-----------------------------------------------------------------------------
     2625
     2626ClpSimplexProgress::~ClpSimplexProgress ()
     2627{
     2628}
     2629// Copy constructor.
     2630ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs)
     2631{
     2632  int i;
     2633  for (i=0;i<CLP_PROGRESS;i++) {
     2634    objective_[i] = rhs.objective_[i];
     2635    infeasibility_[i] = rhs.infeasibility_[i];
     2636    realInfeasibility_[i] = rhs.realInfeasibility_[i];
     2637    numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i];
     2638    iterationNumber_[i]=rhs.iterationNumber_[i];
     2639  }
     2640#ifdef CLP_PROGRESS_WEIGHT
     2641  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
     2642    objectiveWeight_[i] = rhs.objectiveWeight_[i];
     2643    infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
     2644    realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
     2645    numberInfeasibilitiesWeight_[i]=rhs.numberInfeasibilitiesWeight_[i];
     2646    iterationNumberWeight_[i]=rhs.iterationNumberWeight_[i];
     2647  }
     2648  drop_ = rhs.drop_;
     2649  best_ = rhs.best_;
     2650#endif
     2651  initialWeight_ = rhs.initialWeight_;
     2652  for (i=0;i<CLP_CYCLE;i++) {
     2653    //obj_[i]=rhs.obj_[i];
     2654    in_[i]=rhs.in_[i];
     2655    out_[i]=rhs.out_[i];
     2656    way_[i]=rhs.way_[i];
     2657  }
     2658  numberTimes_ = rhs.numberTimes_;
     2659  numberBadTimes_ = rhs.numberBadTimes_;
     2660  model_ = rhs.model_;
     2661  oddState_=rhs.oddState_;
     2662}
     2663// Copy constructor.from model
     2664ClpSimplexProgress::ClpSimplexProgress(ClpSimplex * model)
     2665{
     2666  model_ = model;
     2667  reset();
     2668  initialWeight_=0.0;
     2669}
     2670// Fill from model
     2671void
     2672ClpSimplexProgress::fillFromModel ( ClpSimplex * model )
     2673{
     2674  model_ = model;
     2675  reset();
     2676  initialWeight_=0.0;
     2677}
     2678// Assignment operator. This copies the data
     2679ClpSimplexProgress &
     2680ClpSimplexProgress::operator=(const ClpSimplexProgress & rhs)
     2681{
     2682  if (this != &rhs) {
     2683    int i;
     2684    for (i=0;i<CLP_PROGRESS;i++) {
     2685      objective_[i] = rhs.objective_[i];
     2686      infeasibility_[i] = rhs.infeasibility_[i];
     2687      realInfeasibility_[i] = rhs.realInfeasibility_[i];
     2688      numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i];
     2689      iterationNumber_[i]=rhs.iterationNumber_[i];
     2690    }
     2691#ifdef CLP_PROGRESS_WEIGHT
     2692    for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
     2693      objectiveWeight_[i] = rhs.objectiveWeight_[i];
     2694      infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
     2695      realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
     2696      numberInfeasibilitiesWeight_[i]=rhs.numberInfeasibilitiesWeight_[i];
     2697      iterationNumberWeight_[i]=rhs.iterationNumberWeight_[i];
     2698    }
     2699    drop_ = rhs.drop_;
     2700    best_ = rhs.best_;
     2701#endif
     2702    initialWeight_ = rhs.initialWeight_;
     2703    for (i=0;i<CLP_CYCLE;i++) {
     2704      //obj_[i]=rhs.obj_[i];
     2705      in_[i]=rhs.in_[i];
     2706      out_[i]=rhs.out_[i];
     2707      way_[i]=rhs.way_[i];
     2708    }
     2709    numberTimes_ = rhs.numberTimes_;
     2710    numberBadTimes_ = rhs.numberBadTimes_;
     2711    model_ = rhs.model_;
     2712    oddState_=rhs.oddState_;
     2713  }
     2714  return *this;
     2715}
     2716// Seems to be something odd about exact comparison of doubles on linux
     2717static bool equalDouble(double value1, double value2)
     2718{
     2719
     2720  union { double d; int i[2]; } v1,v2;
     2721  v1.d = value1;
     2722  v2.d = value2;
     2723  if (sizeof(int)*2==sizeof(double))
     2724    return (v1.i[0]==v2.i[0]&&v1.i[1]==v2.i[1]);
     2725  else
     2726    return (v1.i[0]==v2.i[0]);
     2727}
     2728int
     2729ClpSimplexProgress::looping()
     2730{
     2731  if (!model_)
     2732    return -1;
     2733  double objective = model_->rawObjectiveValue();
     2734  double infeasibility;
     2735  double realInfeasibility=0.0;
     2736  int numberInfeasibilities;
     2737  int iterationNumber = model_->numberIterations();
     2738  if (model_->algorithm()<0) {
     2739    // dual
     2740    infeasibility = model_->sumPrimalInfeasibilities();
     2741    numberInfeasibilities = model_->numberPrimalInfeasibilities();
     2742  } else {
     2743    //primal
     2744    infeasibility = model_->sumDualInfeasibilities();
     2745    realInfeasibility = model_->nonLinearCost()->sumInfeasibilities();
     2746    numberInfeasibilities = model_->numberDualInfeasibilities();
     2747  }
     2748  int i;
     2749  int numberMatched=0;
     2750  int matched=0;
     2751  int nsame=0;
     2752  for (i=0;i<CLP_PROGRESS;i++) {
     2753    bool matchedOnObjective = equalDouble(objective,objective_[i]);
     2754    bool matchedOnInfeasibility = equalDouble(infeasibility,infeasibility_[i]);
     2755    bool matchedOnInfeasibilities =
     2756      (numberInfeasibilities==numberInfeasibilities_[i]);
     2757   
     2758    if (matchedOnObjective&&matchedOnInfeasibility&&matchedOnInfeasibilities) {
     2759      matched |= (1<<i);
     2760      // Check not same iteration
     2761      if (iterationNumber!=iterationNumber_[i]) {
     2762        numberMatched++;
     2763        // here mainly to get over compiler bug?
     2764        if (model_->messageHandler()->logLevel()>10)
     2765          printf("%d %d %d %d %d loop check\n",i,numberMatched,
     2766                 matchedOnObjective, matchedOnInfeasibility,
     2767                 matchedOnInfeasibilities);
     2768      } else {
     2769        // stuck but code should notice
     2770        nsame++;
     2771      }
     2772    }
     2773    if (i) {
     2774      objective_[i-1] = objective_[i];
     2775      infeasibility_[i-1] = infeasibility_[i];
     2776      realInfeasibility_[i-1] = realInfeasibility_[i];
     2777      numberInfeasibilities_[i-1]=numberInfeasibilities_[i];
     2778      iterationNumber_[i-1]=iterationNumber_[i];
     2779    }
     2780  }
     2781  objective_[CLP_PROGRESS-1] = objective;
     2782  infeasibility_[CLP_PROGRESS-1] = infeasibility;
     2783  realInfeasibility_[CLP_PROGRESS-1] = realInfeasibility;
     2784  numberInfeasibilities_[CLP_PROGRESS-1]=numberInfeasibilities;
     2785  iterationNumber_[CLP_PROGRESS-1]=iterationNumber;
     2786  if (nsame==CLP_PROGRESS)
     2787    numberMatched=CLP_PROGRESS; // really stuck
     2788  if (model_->progressFlag())
     2789    numberMatched=0;
     2790  numberTimes_++;
     2791  if (numberTimes_<10)
     2792    numberMatched=0;
     2793  // skip if just last time as may be checking something
     2794  if (matched==(1<<(CLP_PROGRESS-1)))
     2795    numberMatched=0;
     2796  if (numberMatched) {
     2797    model_->messageHandler()->message(CLP_POSSIBLELOOP,model_->messages())
     2798      <<numberMatched
     2799      <<matched
     2800      <<numberTimes_
     2801      <<CoinMessageEol;
     2802    numberBadTimes_++;
     2803    if (numberBadTimes_<10) {
     2804      // make factorize every iteration
     2805      model_->forceFactorization(1);
     2806      if (numberBadTimes_<2) {
     2807        startCheck(); // clear other loop check
     2808        if (model_->algorithm()<0) {
     2809          // dual - change tolerance
     2810          model_->setCurrentDualTolerance(model_->currentDualTolerance()*1.05);
     2811          // if infeasible increase dual bound
     2812          if (model_->dualBound()<1.0e17) {
     2813            model_->setDualBound(model_->dualBound()*1.1);
     2814          }
     2815        } else {
     2816          // primal - change tolerance 
     2817          if (numberBadTimes_>3)
     2818            model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance()*1.05);
     2819          // if infeasible increase infeasibility cost
     2820          if (model_->nonLinearCost()->numberInfeasibilities()&&
     2821              model_->infeasibilityCost()<1.0e17) {
     2822            model_->setInfeasibilityCost(model_->infeasibilityCost()*1.1);
     2823          }
     2824        }
     2825      } else {
     2826        // flag
     2827        int iSequence;
     2828        if (model_->algorithm()<0) {
     2829          // dual
     2830          if (model_->dualBound()>1.0e14)
     2831            model_->setDualBound(1.0e14);
     2832          iSequence=in_[CLP_CYCLE-1];
     2833        } else {
     2834          // primal
     2835          if (model_->infeasibilityCost()>1.0e14)
     2836            model_->setInfeasibilityCost(1.0e14);
     2837          iSequence=out_[CLP_CYCLE-1];
     2838        }
     2839        if (iSequence>=0) {
     2840          char x = model_->isColumn(iSequence) ? 'C' :'R';
     2841          if (model_->messageHandler()->logLevel()>=63)
     2842            model_->messageHandler()->message(CLP_SIMPLEX_FLAG,model_->messages())
     2843              <<x<<model_->sequenceWithin(iSequence)
     2844              <<CoinMessageEol;
     2845          // if Gub then needs to be sequenceIn_
     2846          int save=model_->sequenceIn();
     2847          model_->setSequenceIn(iSequence);
     2848          model_->setFlagged(iSequence);
     2849          model_->setSequenceIn(save);
     2850          //printf("flagging %d from loop\n",iSequence);
     2851          startCheck();
     2852        } else {
     2853          // Give up
     2854          if (model_->messageHandler()->logLevel()>=63)
     2855            printf("***** All flagged?\n");
     2856          return 4;
     2857        }
     2858        // reset
     2859        numberBadTimes_=2;
     2860      }
     2861      return -2;
     2862    } else {
     2863      // look at solution and maybe declare victory
     2864      if (infeasibility<1.0e-4) {
     2865        return 0;
     2866      } else {
     2867        model_->messageHandler()->message(CLP_LOOP,model_->messages())
     2868          <<CoinMessageEol;
     2869#ifndef NDEBUG
     2870        printf("debug loop ClpSimplex A\n");
     2871        abort();
     2872#endif
     2873        return 3;
     2874      }
     2875    }
     2876  }
     2877  return -1;
     2878}
     2879// Resets as much as possible
     2880void
     2881ClpSimplexProgress::reset()
     2882{
     2883  int i;
     2884  for (i=0;i<CLP_PROGRESS;i++) {
     2885    if (model_->algorithm()>=0)
     2886      objective_[i] = COIN_DBL_MAX;
     2887    else
     2888      objective_[i] = -COIN_DBL_MAX;
     2889    infeasibility_[i] = -1.0; // set to an impossible value
     2890    realInfeasibility_[i] = COIN_DBL_MAX;
     2891    numberInfeasibilities_[i]=-1;
     2892    iterationNumber_[i]=-1;
     2893  }
     2894#ifdef CLP_PROGRESS_WEIGHT
     2895  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
     2896    objectiveWeight_[i] = COIN_DBL_MAX;
     2897    infeasibilityWeight_[i] = -1.0; // set to an impossible value
     2898    realInfeasibilityWeight_[i] = COIN_DBL_MAX;
     2899    numberInfeasibilitiesWeight_[i]=-1;
     2900    iterationNumberWeight_[i]=-1;
     2901  }
     2902  drop_ =0.0;
     2903  best_ =0.0;
     2904#endif
     2905  for (i=0;i<CLP_CYCLE;i++) {
     2906    //obj_[i]=COIN_DBL_MAX;
     2907    in_[i]=-1;
     2908    out_[i]=-1;
     2909    way_[i]=0;
     2910  }
     2911  numberTimes_ = 0;
     2912  numberBadTimes_ = 0;
     2913  oddState_=0;
     2914}
     2915// Returns previous objective (if -1) - current if (0)
     2916double
     2917ClpSimplexProgress::lastObjective(int back) const
     2918{
     2919  return objective_[CLP_PROGRESS-1-back];
     2920}
     2921// Returns previous infeasibility (if -1) - current if (0)
     2922double
     2923ClpSimplexProgress::lastInfeasibility(int back) const
     2924{
     2925  return realInfeasibility_[CLP_PROGRESS-1-back];
     2926}
     2927// Sets real primal infeasibility
     2928void
     2929ClpSimplexProgress::setInfeasibility(double value)
     2930{
     2931  for (int i=1;i<CLP_PROGRESS;i++)
     2932    realInfeasibility_[i-1] = realInfeasibility_[i];
     2933  realInfeasibility_[CLP_PROGRESS-1]=value;
     2934}
     2935// Modify objective e.g. if dual infeasible in dual
     2936void
     2937ClpSimplexProgress::modifyObjective(double value)
     2938{
     2939  objective_[CLP_PROGRESS-1]=value;
     2940}
     2941// Returns previous iteration number (if -1) - current if (0)
     2942int
     2943ClpSimplexProgress::lastIterationNumber(int back) const
     2944{
     2945  return iterationNumber_[CLP_PROGRESS-1-back];
     2946}
     2947// clears iteration numbers (to switch off panic)
     2948void
     2949ClpSimplexProgress::clearIterationNumbers()
     2950{
     2951  for (int i=0;i<CLP_PROGRESS;i++)
     2952    iterationNumber_[i]=-1;
     2953}
     2954// Start check at beginning of whileIterating
     2955void
     2956ClpSimplexProgress::startCheck()
     2957{
     2958  int i;
     2959  for (i=0;i<CLP_CYCLE;i++) {
     2960    //obj_[i]=COIN_DBL_MAX;
     2961    in_[i]=-1;
     2962    out_[i]=-1;
     2963    way_[i]=0;
     2964  }
     2965}
     2966// Returns cycle length in whileIterating
     2967int
     2968ClpSimplexProgress::cycle(int in, int out,int wayIn,int wayOut)
     2969{
     2970  int i;
     2971#if 0
     2972  if (model_->numberIterations()>206571) {
     2973    printf("in %d out %d\n",in,out);
     2974    for (i=0;i<CLP_CYCLE;i++)
     2975      printf("cy %d in %d out %d\n",i,in_[i],out_[i]);
     2976  }
     2977#endif
     2978  int matched=0;
     2979  // first see if in matches any out
     2980  for (i=1;i<CLP_CYCLE;i++) {
     2981    if (in==out_[i]) {
     2982      // even if flip then suspicious
     2983      matched=-1;
     2984      break;
     2985    }
     2986  }
     2987#if 0
     2988  if (!matched||in_[0]<0) {
     2989    // can't be cycle
     2990    for (i=0;i<CLP_CYCLE-1;i++) {
     2991      //obj_[i]=obj_[i+1];
     2992      in_[i]=in_[i+1];
     2993      out_[i]=out_[i+1];
     2994      way_[i]=way_[i+1];
     2995    }
     2996  } else {
     2997    // possible cycle
     2998    matched=0;
     2999    for (i=0;i<CLP_CYCLE-1;i++) {
     3000      int k;
     3001      char wayThis = way_[i];
     3002      int inThis = in_[i];
     3003      int outThis = out_[i];
     3004      //double objThis = obj_[i];
     3005      for(k=i+1;k<CLP_CYCLE;k++) {
     3006        if (inThis==in_[k]&&outThis==out_[k]&&wayThis==way_[k]) {
     3007          int distance = k-i;
     3008          if (k+distance<CLP_CYCLE) {
     3009            // See if repeats
     3010            int j=k+distance;
     3011            if (inThis==in_[j]&&outThis==out_[j]&&wayThis==way_[j]) {
     3012              matched=distance;
     3013              break;
     3014            }
     3015          } else {
     3016            matched=distance;
     3017            break;
     3018          }
     3019        }
     3020      }
     3021      //obj_[i]=obj_[i+1];
     3022      in_[i]=in_[i+1];
     3023      out_[i]=out_[i+1];
     3024      way_[i]=way_[i+1];
     3025    }
     3026  }
     3027#else
     3028  if (matched&&in_[0]>=0) {
     3029    // possible cycle - only check [0] against all
     3030    matched=0;
     3031    int nMatched=0;
     3032    char way0 = way_[0];
     3033    int in0 = in_[0];
     3034    int out0 = out_[0];
     3035    //double obj0 = obj_[i];
     3036    for(int k=1;k<CLP_CYCLE-4;k++) {
     3037      if (in0==in_[k]&&out0==out_[k]&&way0==way_[k]) {
     3038        nMatched++;
     3039        // See if repeats
     3040        int end = CLP_CYCLE-k;
     3041        int j;
     3042        for ( j=1;j<end;j++) {
     3043          if (in_[j+k]!=in_[j]||out_[j+k]!=out_[j]||way_[j+k]!=way_[j])
     3044            break;
     3045        }
     3046        if (j==end) {
     3047          matched=k;
     3048          break;
     3049        }
     3050      }
     3051    }
     3052    // If three times then that is too much even if not regular
     3053    if (matched<=0&&nMatched>1)
     3054      matched=100;
     3055  }
     3056  for (i=0;i<CLP_CYCLE-1;i++) {
     3057    //obj_[i]=obj_[i+1];
     3058    in_[i]=in_[i+1];
     3059    out_[i]=out_[i+1];
     3060    way_[i]=way_[i+1];
     3061  }
     3062#endif
     3063  char way = 1-wayIn+4*(1-wayOut);
     3064  //obj_[i]=model_->objectiveValue();
     3065  in_[CLP_CYCLE-1]=in;
     3066  out_[CLP_CYCLE-1]=out;
     3067  way_[CLP_CYCLE-1]=way;
     3068  return matched;
     3069}
  • trunk/Clp/src/ClpSolve.hpp

    r1055 r1147  
    204204  //@}
    205205};
     206
     207/// For saving extra information to see if looping.
     208class ClpSimplexProgress {
     209
     210public:
     211
     212
     213  /**@name Constructors and destructor and copy */
     214  //@{
     215  /// Default constructor
     216    ClpSimplexProgress (  );
     217
     218  /// Constructor from model
     219    ClpSimplexProgress ( ClpSimplex * model );
     220
     221  /// Copy constructor.
     222  ClpSimplexProgress(const ClpSimplexProgress &);
     223
     224  /// Assignment operator. This copies the data
     225    ClpSimplexProgress & operator=(const ClpSimplexProgress & rhs);
     226  /// Destructor
     227   ~ClpSimplexProgress (  );
     228  /// Resets as much as possible
     229   void reset();
     230  /// Fill from model
     231    void fillFromModel ( ClpSimplex * model );
     232
     233  //@}
     234
     235  /**@name Check progress */
     236  //@{
     237  /** Returns -1 if okay, -n+1 (n number of times bad) if bad but action taken,
     238      >=0 if give up and use as problem status
     239  */
     240    int looping (  );
     241  /// Start check at beginning of whileIterating
     242  void startCheck();
     243  /// Returns cycle length in whileIterating
     244  int cycle(int in, int out,int wayIn,int wayOut);
     245
     246  /// Returns previous objective (if -1) - current if (0)
     247  double lastObjective(int back=1) const;
     248  /// Set real primal infeasibility and move back
     249  void setInfeasibility(double value);
     250  /// Returns real primal infeasibility (if -1) - current if (0)
     251  double lastInfeasibility(int back=1) const;
     252  /// Modify objective e.g. if dual infeasible in dual
     253  void modifyObjective(double value);
     254  /// Returns previous iteration number (if -1) - current if (0)
     255  int lastIterationNumber(int back=1) const;
     256  /// clears all iteration numbers (to switch off panic)
     257  void clearIterationNumbers();
     258  /// Odd state
     259  inline void newOddState()
     260  { oddState_= - oddState_-1;}
     261  inline void endOddState()
     262  { oddState_=abs(oddState_);}
     263  inline void clearOddState()
     264  { oddState_=0;}
     265  inline int oddState() const
     266  { return oddState_;}
     267  /// number of bad times
     268  inline int badTimes() const
     269  { return numberBadTimes_;}
     270  inline void clearBadTimes()
     271  { numberBadTimes_=0;}
     272
     273  //@}
     274  /**@name Data  */
     275#define CLP_PROGRESS 5
     276  //#define CLP_PROGRESS_WEIGHT 10
     277  //@{
     278  /// Objective values
     279  double objective_[CLP_PROGRESS];
     280  /// Sum of infeasibilities for algorithm
     281  double infeasibility_[CLP_PROGRESS];
     282  /// Sum of real primal infeasibilities for primal
     283  double realInfeasibility_[CLP_PROGRESS];
     284#ifdef CLP_PROGRESS_WEIGHT
     285  /// Objective values for weights
     286  double objectiveWeight_[CLP_PROGRESS_WEIGHT];
     287  /// Sum of infeasibilities for algorithm for weights
     288  double infeasibilityWeight_[CLP_PROGRESS_WEIGHT];
     289  /// Sum of real primal infeasibilities for primal for weights
     290  double realInfeasibilityWeight_[CLP_PROGRESS_WEIGHT];
     291  /// Drop  for weights
     292  double drop_;
     293  /// Best? for weights
     294  double best_;
    206295#endif
     296  /// Initial weight for weights
     297  double initialWeight_;
     298#define CLP_CYCLE 12
     299  /// For cycle checking
     300  //double obj_[CLP_CYCLE];
     301  int in_[CLP_CYCLE];
     302  int out_[CLP_CYCLE];
     303  char way_[CLP_CYCLE];
     304  /// Pointer back to model so we can get information
     305  ClpSimplex * model_;
     306  /// Number of infeasibilities
     307  int numberInfeasibilities_[CLP_PROGRESS];
     308  /// Iteration number at which occurred
     309  int iterationNumber_[CLP_PROGRESS];
     310#ifdef CLP_PROGRESS_WEIGHT
     311  /// Number of infeasibilities for weights
     312  int numberInfeasibilitiesWeight_[CLP_PROGRESS_WEIGHT];
     313  /// Iteration number at which occurred for weights
     314  int iterationNumberWeight_[CLP_PROGRESS_WEIGHT];
     315#endif
     316  /// Number of times checked (so won't stop too early)
     317  int numberTimes_;
     318  /// Number of times it looked like loop
     319  int numberBadTimes_;
     320  /// If things are in an odd state
     321  int oddState_;
     322  //@}
     323};
     324#endif
Note: See TracChangeset for help on using the changeset viewer.