Changeset 208 for branches


Ignore:
Timestamp:
Sep 12, 2003 5:20:25 PM (16 years ago)
Author:
forrest
Message:

Trying to make faster

Location:
branches/pre
Files:
1 added
18 edited

Legend:

Unmodified
Added
Removed
  • branches/pre/ClpDualRowSteepest.cpp

    r205 r208  
    212212  }
    213213  int numberWanted;
    214   if (mode_<2)
     214  if (mode_<2 ) {
    215215    numberWanted = number+1;
    216   else
     216  } else if (mode_==2) {
    217217    numberWanted = max(2000,number/8);
     218  } else {
     219    int numberElements = model_->factorization()->numberElements();
     220    double ratio = (double) numberElements/(double) model_->numberRows();
     221    numberWanted = max(2000,number/8);
     222    if (ratio<1.0) {
     223      numberWanted = max(2000,number/20);
     224    } else if (ratio>10.0) {
     225      ratio = number * (ratio/80.0);
     226      if (ratio>number)
     227        numberWanted=number+1;
     228      else
     229        numberWanted = max(2000,(int) ratio);
     230    }
     231  }
    218232  int iPass;
    219233  // Setup two passes
     
    225239  start[3]=start[0];
    226240  //double largestWeight=0.0;
     241  //double smallestWeight=1.0e100;
    227242  for (iPass=0;iPass<2;iPass++) {
    228243    int end = start[2*iPass+1];
     
    231246      double value = infeas[iRow];
    232247      if (value>tolerance) {
     248        //#define OUT_EQ
     249#ifdef OUT_EQ
     250        {
     251          int iSequence = pivotVariable[iRow];
     252          if (upper[iSequence]==lower[iSequence])
     253            value *= 2.0;
     254        }
     255#endif
    233256        double weight = weights_[iRow];
     257        //largestWeight = max(largestWeight,weight);
     258        //smallestWeight = min(smallestWeight,weight);
    234259        //double dubious = dubiousWeights_[iRow];
    235260        //weight *= dubious;
     
    270295      break;
    271296  }
     297  //printf("smallest %g largest %g\n",smallestWeight,largestWeight);
    272298  return chosenRow;
    273299}
  • branches/pre/ClpMessage.cpp

    r190 r208  
    2222  {CLP_SIMPLEX_ACCURACY,6,3,"Primal error %g, dual error %g"},
    2323  {CLP_SIMPLEX_BADFACTOR,7,1,"Singular factorization of basis - status %d"},
    24   {CLP_SIMPLEX_BOUNDTIGHTEN,8,1,"Bounds were tightened %d times"},
     24  {CLP_SIMPLEX_BOUNDTIGHTEN,8,3,"Bounds were tightened %d times"},
    2525  {CLP_SIMPLEX_INFEASIBILITIES,9,1,"%d infeasibilities"},
    2626  {CLP_SIMPLEX_FLAG,10,3,"Flagging variable %c%d"},
     
    2828  {CLP_DUAL_CHECKB,12,2,"New dual bound of %g"},
    2929  {CLP_DUAL_ORIGINAL,13,3,"Going back to original objective"},
    30   {CLP_SIMPLEX_PERTURB,14,1,"Perturbing problem by %g"},
     30  {CLP_SIMPLEX_PERTURB,14,1,"Perturbing problem by %g %% of %g - largest change %g (%% %g) - largest zero change %g"},
    3131  {CLP_PRIMAL_ORIGINAL,15,2,"Going back to original tolerance"},
    3232  {CLP_PRIMAL_WEIGHT,16,2,"New infeasibility weight of %g"},
     
    6363  {CLP_QUADRATIC_BOTH,108,32,"%s %d (%g) and %d (%g) both basic"},
    6464  {CLP_QUADRATIC_PRIMAL_DETAILS,109,32,"coeff %g, %g, %g - dj %g - deriv zero at %g, sj at %g"},
     65  {CLP_IDIOT_ITERATION,30,1,"%d infeas %g, obj %g - mu %g, its %d, %d interior"},
     66  {CLP_INFEASIBLE,3003,1,"Analysis indicates model infeasible"},
     67  {CLP_MATRIX_CHANGE,31,1,"Matrix can not be converted into %s"},
     68  {CLP_TIMING,32,1,"%s objective %.9g - %d iterations time %.2f2%?, Presolve %.2f%?, Idiot %.2f%?"},
     69  {CLP_INTERVAL_TIMING,33,2,"%s took %.2f seconds (total %.2f)"},
    6570  {CLP_DUMMY_END,999999,0,""}
    6671};
  • branches/pre/ClpModel.cpp

    r205 r208  
    6666  strParam_[ClpProbName] = "ClpDefaultName";
    6767  handler_ = new CoinMessageHandler();
    68   handler_->setLogLevel(2);
     68  handler_->setLogLevel(1);
    6969  messages_ = ClpMessage();
    7070}
  • branches/pre/ClpNetworkMatrix.cpp

    r207 r208  
    370370    dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
    371371  bool packed = rowArray->packedMode();
    372   if (numberInRowArray>0.3*numberRows||!rowCopy) {
     372  double factor = 0.3;
     373  // We may not want to do by row if there may be cache problems
     374  int numberColumns = model->numberColumns();
     375  // It would be nice to find L2 cache size - for moment 512K
     376  // Be slightly optimistic
     377  if (numberColumns*sizeof(double)>1000000) {
     378    if (numberRows*10<numberColumns)
     379      factor=0.1;
     380    else if (numberRows*4<numberColumns)
     381      factor=0.15;
     382    else if (numberRows*2<numberColumns)
     383      factor=0.2;
     384    if (model->numberIterations()%50==0)
     385      printf("%d nonzero\n",numberInRowArray);
     386  }
     387  if (numberInRowArray>factor*numberRows||!rowCopy) {
    373388    // do by column
    374389    int iColumn;
  • branches/pre/ClpPackedMatrix.cpp

    r205 r208  
    270270    dynamic_cast< ClpPackedMatrix*>(model->rowCopy());
    271271  bool packed = rowArray->packedMode();
    272   if (numberInRowArray>0.333*numberRows||!rowCopy) {
     272  double factor = 0.3;
     273  // We may not want to do by row if there may be cache problems
     274  int numberColumns = model->numberColumns();
     275  // It would be nice to find L2 cache size - for moment 512K
     276  // Be slightly optimistic
     277  if (numberColumns*sizeof(double)>1000000) {
     278    if (numberRows*10<numberColumns)
     279      factor=0.1;
     280    else if (numberRows*4<numberColumns)
     281      factor=0.15;
     282    else if (numberRows*2<numberColumns)
     283      factor=0.2;
     284    if (model->numberIterations()%50==0)
     285      printf("%d nonzero\n",numberInRowArray);
     286  }
     287  if (numberInRowArray>factor*numberRows||!rowCopy) {
    273288    // do by column
    274289    int iColumn;
     
    11661181    ClpFillN ( rowScale, numberRows,1.0);
    11671182    ClpFillN ( columnScale, numberColumns,1.0);
    1168     int numberPass=3;
    1169     double overallLargest;
    1170     double overallSmallest;
    1171     while (numberPass) {
    1172       overallLargest=0.0;
    1173       overallSmallest=1.0e50;
    1174       numberPass--;
    1175       // Geometric mean on row scales
     1183    double overallLargest=-1.0e-30;
     1184    double overallSmallest=1.0e30;
     1185    if (model->scalingFlag()==1) {
     1186      // Maximum in each row
    11761187      for (iRow=0;iRow<numberRows;iRow++) {
    11771188        if (usefulRow[iRow]) {
    11781189          CoinBigIndex j;
    11791190          largest=1.0e-20;
    1180           smallest=1.0e50;
    11811191          for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
    11821192            int iColumn = column[j];
    11831193            if (usefulColumn[iColumn]) {
    11841194              double value = fabs(element[j]);
     1195              largest = max(largest,value);
     1196            }
     1197          }
     1198          rowScale[iRow]=1.0/largest;
     1199          overallLargest = max(overallLargest,largest);
     1200          overallSmallest = min(overallSmallest,largest);
     1201        }
     1202      }
     1203    } else {
     1204      assert(model->scalingFlag()==2);
     1205      int numberPass=3;
     1206#ifdef USE_OBJECTIVE
     1207      // This will be used to help get scale factors
     1208      double * objective = new double[numberColumns];
     1209      memcpy(objective,model->costRegion(1),numberColumns*sizeof(double));
     1210      double objScale=1.0;
     1211#endif
     1212      while (numberPass) {
     1213        overallLargest=0.0;
     1214        overallSmallest=1.0e50;
     1215        numberPass--;
     1216        // Geometric mean on row scales
     1217        for (iRow=0;iRow<numberRows;iRow++) {
     1218          if (usefulRow[iRow]) {
     1219            CoinBigIndex j;
     1220            largest=1.0e-20;
     1221            smallest=1.0e50;
     1222            for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     1223              int iColumn = column[j];
     1224              if (usefulColumn[iColumn]) {
     1225                double value = fabs(element[j]);
     1226                // Don't bother with tiny elements
     1227                if (value>1.0e-30) {
     1228                  value *= columnScale[iColumn];
     1229                  largest = max(largest,value);
     1230                  smallest = min(smallest,value);
     1231                }
     1232              }
     1233            }
     1234            rowScale[iRow]=1.0/sqrt(smallest*largest);
     1235            overallLargest = max(largest*rowScale[iRow],overallLargest);
     1236            overallSmallest = min(smallest*rowScale[iRow],overallSmallest);
     1237          }
     1238        }
     1239#ifdef USE_OBJECTIVE
     1240        largest=1.0e-20;
     1241        smallest=1.0e50;
     1242        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1243          if (usefulColumn[iColumn]) {
     1244            double value = fabs(objective[iColumn]);
     1245            // Don't bother with tiny elements
     1246            if (value>1.0e-30) {
     1247              value *= columnScale[iColumn];
     1248              largest = max(largest,value);
     1249              smallest = min(smallest,value);
     1250            }
     1251          }
     1252        }
     1253        objScale=1.0/sqrt(smallest*largest);
     1254#endif
     1255        model->messageHandler()->message(CLP_PACKEDSCALE_WHILE,*model->messagesPointer())
     1256          <<overallSmallest
     1257          <<overallLargest
     1258          <<CoinMessageEol;
     1259        // skip last column round
     1260        if (numberPass==1)
     1261          break;
     1262        // Geometric mean on column scales
     1263        for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1264          if (usefulColumn[iColumn]) {
     1265            CoinBigIndex j;
     1266            largest=1.0e-20;
     1267            smallest=1.0e50;
     1268            for (j=columnStart[iColumn];
     1269                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1270              iRow=row[j];
     1271              double value = fabs(elementByColumn[j]);
    11851272              // Don't bother with tiny elements
    1186               if (value>1.0e-30) {
    1187                 value *= columnScale[iColumn];
     1273              if (value>1.0e-30&&usefulRow[iRow]) {
     1274                value *= rowScale[iRow];
    11881275                largest = max(largest,value);
    11891276                smallest = min(smallest,value);
    11901277              }
    11911278            }
    1192           }
    1193           rowScale[iRow]=1.0/sqrt(smallest*largest);
    1194           overallLargest = max(largest*rowScale[iRow],overallLargest);
    1195           overallSmallest = min(smallest*rowScale[iRow],overallSmallest);
    1196         }
    1197       }
    1198       model->messageHandler()->message(CLP_PACKEDSCALE_WHILE,*model->messagesPointer())
    1199         <<overallSmallest
    1200         <<overallLargest
    1201         <<CoinMessageEol;
    1202       // skip last column round
    1203       if (numberPass==1)
    1204         break;
    1205       // Geometric mean on column scales
    1206       for (iColumn=0;iColumn<numberColumns;iColumn++) {
    1207         if (usefulColumn[iColumn]) {
    1208           CoinBigIndex j;
    1209           largest=1.0e-20;
    1210           smallest=1.0e50;
    1211           for (j=columnStart[iColumn];
    1212                j<columnStart[iColumn]+columnLength[iColumn];j++) {
    1213             iRow=row[j];
    1214             double value = fabs(elementByColumn[j]);
    1215             // Don't bother with tiny elements
    1216             if (value>1.0e-30&&usefulRow[iRow]) {
    1217               value *= rowScale[iRow];
     1279#ifdef USE_OBJECTIVE
     1280            if (fabs(objective[iColumn])>1.0e-30) {
     1281              double value = fabs(objective[iColumn])*objScale;
    12181282              largest = max(largest,value);
    12191283              smallest = min(smallest,value);
    12201284            }
    1221           }
    1222           columnScale[iColumn]=1.0/sqrt(smallest*largest);
    1223         }
    1224       }
    1225     }
    1226     // final pass to scale columns so largest is 1.0
    1227     overallLargest=0.0;
     1285#endif
     1286            columnScale[iColumn]=1.0/sqrt(smallest*largest);
     1287          }
     1288        }
     1289      }
     1290#ifdef USE_OBJECTIVE
     1291      delete [] objective;
     1292      printf("obj scale %g - use it if you want\n",objScale);
     1293#endif
     1294    }
     1295    // final pass to scale columns so largest is reasonable
     1296    // See what smallest will be if largest is 1.0
    12281297    overallSmallest=1.0e50;
    12291298    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     
    12411310          }
    12421311        }
    1243         columnScale[iColumn]=1.0/largest;
    1244         overallLargest = max(overallLargest,largest*columnScale[iColumn]);
     1312        if (overallSmallest*largest>smallest)
     1313          overallSmallest = smallest/largest;
     1314      }
     1315    }
     1316    overallLargest=1.0;
     1317    if (overallSmallest<1.0e-3)
     1318      overallLargest = 1.0/sqrt(overallSmallest);
     1319    overallLargest = min(1000.0,overallLargest);
     1320    overallSmallest=1.0e50;
     1321    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1322      if (usefulColumn[iColumn]) {
     1323        CoinBigIndex j;
     1324        largest=1.0e-20;
     1325        smallest=1.0e50;
     1326        for (j=columnStart[iColumn];
     1327             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1328          iRow=row[j];
     1329          if(elementByColumn[j]&&usefulRow[iRow]) {
     1330            double value = fabs(elementByColumn[j]*rowScale[iRow]);
     1331            largest = max(largest,value);
     1332            smallest = min(smallest,value);
     1333          }
     1334        }
     1335        columnScale[iColumn]=overallLargest/largest;
    12451336        overallSmallest = min(overallSmallest,smallest*columnScale[iColumn]);
    12461337      }
     
    12501341      <<overallLargest
    12511342      <<CoinMessageEol;
    1252    
    12531343    delete [] usefulRow;
    12541344    delete [] usefulColumn;
  • branches/pre/ClpPlusMinusOneMatrix.cpp

    r207 r208  
    459459  ClpPlusMinusOneMatrix* rowCopy =
    460460    dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
    461   if (numberInRowArray>0.3*numberRows||!rowCopy) {
     461  double factor = 0.3;
     462  // We may not want to do by row if there may be cache problems
     463  int numberColumns = model->numberColumns();
     464  // It would be nice to find L2 cache size - for moment 512K
     465  // Be slightly optimistic
     466  if (numberColumns*sizeof(double)>1000000) {
     467    if (numberRows*10<numberColumns)
     468      factor=0.1;
     469    else if (numberRows*4<numberColumns)
     470      factor=0.15;
     471    else if (numberRows*2<numberColumns)
     472      factor=0.2;
     473    if (model->numberIterations()%50==0)
     474      printf("%d nonzero\n",numberInRowArray);
     475  }
     476  if (numberInRowArray>factor*numberRows||!rowCopy) {
    462477    assert (!y->getNumElements());
    463478    // do by column
  • branches/pre/ClpSimplex.cpp

    r205 r208  
    19291929  bool goodMatrix=true;
    19301930  int saveLevel=handler_->logLevel();
    1931   if (problemStatus_==10)
     1931  if (problemStatus_==10) {
    19321932    handler_->setLogLevel(0); // switch off messages
     1933  } else if (factorization_) {
     1934    // match up factorization messages
     1935    if (handler_->logLevel()<3)
     1936      factorization_->messageLevel(0);
     1937    else
     1938      factorization_->messageLevel(3);
     1939  }
    19331940  int i;
    19341941  if ((what&(16+32))!=0) {
     
    22872294ClpSimplex::setPerturbation(int value)
    22882295{
    2289   if(value<=100&&value >=-50) {
     2296  if(value<=100&&value >=-1000) {
    22902297    perturbation_=value;
    22912298  }
     
    41864193    int saveThreshold = factorization_->denseThreshold();
    41874194    factorization_->setDenseThreshold(0);
    4188 
     4195#if 0
    41894196    // do perturbation if asked for
    41904197
     
    41964203      }
    41974204    }
     4205#endif
    41984206    // for primal we will change bounds using infeasibilityCost_
    41994207    if (nonLinearCost_==NULL&&algorithm_>0) {
  • branches/pre/ClpSimplexDual.cpp

    r205 r208  
    221221    // If values pass then scale pi
    222222    if (ifValuesPass) {
     223      if (!problemStatus_&&perturbation_<100)
     224        perturb();
    223225      int i;
    224226      if (scalingFlag_>0) {
     
    250252        }
    251253      }
    252     }
     254    } 
    253255
    254256    double objectiveChange;
     
    258260    int lastCleaned=0; // last time objective or bounds cleaned up
    259261
    260    
     262    if (!ifValuesPass) {
     263      // Check optimal
     264      if (!numberDualInfeasibilities_&&!numberPrimalInfeasibilities_)
     265        problemStatus_=0;
     266    }
     267    if (problemStatus_<0&&perturbation_<100) {
     268      perturb();
     269      // Can't get here if values pass
     270      gutsOfSolution(NULL,NULL);
     271      if (handler_->logLevel()>2) {
     272        handler_->message(CLP_SIMPLEX_STATUS,messages_)
     273          <<numberIterations_<<objectiveValue();
     274        handler_->printing(sumPrimalInfeasibilities_>0.0)
     275          <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
     276        handler_->printing(sumDualInfeasibilities_>0.0)
     277          <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
     278        handler_->printing(numberDualInfeasibilitiesWithoutFree_
     279                           <numberDualInfeasibilities_)
     280                             <<numberDualInfeasibilitiesWithoutFree_;
     281        handler_->message()<<CoinMessageEol;
     282      }
     283    }
     284       
    261285    // This says whether to restore things etc
    262286    int factorType=0;
     
    19611985              bool take=false;
    19621986              double absAlpha = fabs(alpha);
     1987              // User could do anything to break ties here
    19631988              double weight;
    19641989              if (dubiousWeights)
     
    29062931  double perturbation=1.0e-20;
    29072932  // maximum fraction of cost to perturb
    2908   double maximumFraction = 1.0e-4;
    2909   double factor=1.0e-8;
     2933  double maximumFraction = 1.0e-5;
     2934  double constantPerturbation = 100.0*dualTolerance_;
     2935  const int * lengths = matrix_->getVectorLengths();
     2936  int maxLength=0;
     2937  int minLength=numberRows_;
     2938  if (perturbation_==50) {
     2939    int count[101];
     2940    int numberTotal=0;
     2941    memset(count,0,101*sizeof(int));
     2942    int iColumn;
     2943    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     2944      if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
     2945        double value = fabs(objectiveWork_[iColumn]);
     2946        int k;
     2947        numberTotal++;
     2948        if (!value) {
     2949          k=0;
     2950        } else if (value>99.5) {
     2951          k=100;
     2952        } else {
     2953          double lo = floor(value+0.5);
     2954          if (fabs(lo-value)<1.0e-8)
     2955            k = (int) lo;
     2956          else
     2957            k=100;
     2958        }
     2959        count[k]++;
     2960        int length = lengths[iColumn];
     2961        if (length>2) {
     2962          maxLength = max(maxLength,length);
     2963          minLength = min(minLength,length);
     2964        }
     2965      }
     2966    }
     2967    int numberZero = count[0];
     2968    int numberNonZero = numberTotal-numberZero;
     2969    int numberOdd=count[100];
     2970    int numberInteger = numberNonZero-numberOdd;
     2971    if (3*numberOdd>numberTotal) {
     2972      // be more subtle
     2973      maximumFraction=1.0e-6;
     2974      constantPerturbation *= 0.1;
     2975    } else if (count[1]==numberInteger&&!numberOdd) {
     2976      // be less subtle
     2977      maximumFraction=2.0e-5;
     2978    }
     2979  } else {
     2980    int iColumn;
     2981    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     2982      if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
     2983        int length = lengths[iColumn];
     2984        if (length>2) {
     2985          maxLength = max(maxLength,length);
     2986          minLength = min(minLength,length);
     2987        }
     2988      }
     2989    }
     2990  }
    29102991  // If > 70 then do rows
    2911   if (perturbation_>70) {
     2992  if (perturbation_>=70) {
    29122993    modifyRowCosts=true;
    29132994    perturbation_ -= 20;
    2914     //printf("Row costs modified, ");
    2915   }
     2995    printf("Row costs modified, ");
     2996  }
     2997  bool uniformChange=false;
    29162998  if (perturbation_>50) {
    29172999    // Experiment
    2918     // factor could be 1.0e-8 to 1.0 (51,56,61,66 are 1.0e-8)
    2919     // maximumFraction could be 1.0e-6 to 1.0e-3 (51-55 are 1.0e-6)
    2920     double f[]={1.0e-8,1.0e-6,1.0e-4,1.0e-2,1.0};
    2921     double m[]={1.0e-6,1.0e-5,1.0e-4,1.0e-3};
    2922     factor = f[((perturbation_-1)%5)];
    2923     maximumFraction = m[(perturbation_-51)/5];
     3000    // maximumFraction could be 1.0e-10 to 1.0
     3001    double m[]={1.0e-10,1.0e-9,1.0e-8,1.0e-7,1.0e-6,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,1.0};
     3002    maximumFraction = m[min(perturbation_-51,10)];
    29243003  }
    29253004  int iRow,iColumn;
     3005  double smallestNonZero=1.0e100;
    29263006  if (perturbation_>=50) {
    29273007    perturbation = 1.0e-8;
     
    29303010        double value = fabs(rowObjectiveWork_[iRow]);
    29313011        perturbation = max(perturbation,value);
    2932         if (value)
     3012        if (value) {
    29333013          modifyRowCosts=true;
     3014          smallestNonZero = min(smallestNonZero,value);
     3015        }
    29343016      }
    29353017    }
     
    29393021          fabs(objectiveWork_[iColumn]);
    29403022        perturbation = max(perturbation,value);
    2941       }
    2942     }
    2943     perturbation *= factor;
     3023        if (value) {
     3024          smallestNonZero = min(smallestNonZero,value);
     3025        }
     3026      }
     3027    }
     3028    perturbation = min(perturbation,smallestNonZero/maximumFraction);
    29443029  } else {
    29453030    // user is in charge
    2946     maximumFraction = 1.0e100;
     3031    maximumFraction = 1.0e-1;
    29473032    // but some experiments
    29483033    if (perturbation_<=-900) {
    29493034      modifyRowCosts=true;
    29503035      perturbation_ += 1000;
    2951       //printf("Row costs modified, ");
    2952     }
    2953     if (perturbation_<-50) {
     3036      printf("Row costs modified, ");
     3037    }
     3038    if (perturbation_<=-10) {
     3039      perturbation_ += 10;
    29543040      maximumFraction = 1.0;
    2955       while (perturbation_<-50) {
     3041      if ((-perturbation_)%100>=10) {
     3042        uniformChange=true;
     3043        perturbation_+=20;
     3044      }
     3045      while (perturbation_<-10) {
    29563046        perturbation_ += 100;
    29573047        maximumFraction *= 1.0e-1;
     
    29603050    perturbation = pow(10.0,perturbation_);
    29613051  }
    2962   //printf("factor %g, maximumFraction %g\n",factor,maximumFraction);
     3052  double largestZero=0.0;
     3053  double largest=0.0;
     3054  double largestPerCent=0.0;
    29633055  // modify costs
    2964   handler_->message(CLP_SIMPLEX_PERTURB,messages_)
    2965     <<perturbation
    2966     <<CoinMessageEol;
     3056  bool printOut=(handler_->logLevel()==63);
    29673057  if (modifyRowCosts) {
    29683058    for (iRow=0;iRow<numberRows_;iRow++) {
     
    29703060        double value = perturbation;
    29713061        double currentValue = rowObjectiveWork_[iRow];
    2972         value = min(value,maximumFraction*fabs(currentValue)+1.0e-6);
     3062        value = min(value,maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-3));
    29733063        if (rowLowerWork_[iRow]>-largeValue_) {
    29743064          if (fabs(rowLowerWork_[iRow])<fabs(rowUpperWork_[iRow]))
     
    29813071          value=0.0;
    29823072        }
     3073        if (currentValue) {
     3074          largest = max(largest,fabs(value));
     3075          if (fabs(value)>fabs(currentValue)*largestPerCent)
     3076            largestPerCent=fabs(value/currentValue);
     3077        } else {
     3078          largestZero=max(largestZero,fabs(value));
     3079        }
     3080        if (printOut)
     3081          printf("row %d cost %g change %g\n",iRow,rowObjectiveWork_[iRow],value);
    29833082        rowObjectiveWork_[iRow] += value;
    29843083      }
    29853084    }
    29863085  }
     3086  // more its but faster double weight[]={1.0e-4,1.0e-2,1.0e-1,1.0,2.0,10.0,100.0,200.0,400.0,600.0,1000.0};
     3087  // good its double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,5.0,10.0,20.0,30.0,40.0,100.0};
     3088  double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,5.0,10.0,20.0,30.0,40.0,100.0};
     3089  //double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,20.0,50.0,100.0,120.0,130.0,140.0,200.0};
     3090  double extraWeight=10.0;
     3091  // Scale back if wanted
     3092  double weight2[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0};
     3093  if (constantPerturbation<99.0*dualTolerance_) {
     3094    perturbation *= 0.1;
     3095    extraWeight=0.5;
     3096    memcpy(weight,weight2,sizeof(weight2));
     3097  }
     3098  // adjust weights if all columns long
     3099  double factor=1.0;
     3100  if (maxLength) {
     3101    factor = 3.0/(double) minLength;
     3102  }
     3103  // Make variables with more elements more expensive
     3104  const double m1 = 0.5;
    29873105  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    29883106    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
    29893107      double value = perturbation;
    29903108      double currentValue = objectiveWork_[iColumn];
    2991       value = min(value,maximumFraction*fabs(currentValue)+1.0e-6);
     3109      value = min(value,constantPerturbation+maximumFraction*(fabs(currentValue)+1.0e-1*perturbation+1.0e-8));
     3110      //value = min(value,constantPerturbation;+maximumFraction*fabs(currentValue));
     3111      double value2 = constantPerturbation+1.0e-1*smallestNonZero;
     3112      if (uniformChange) {
     3113        value = maximumFraction;
     3114        value2=maximumFraction;
     3115      }
    29923116      if (columnLowerWork_[iColumn]>-largeValue_) {
    29933117        if (fabs(columnLowerWork_[iColumn])<
    2994             fabs(columnUpperWork_[iColumn]))
    2995           value *= CoinDrand48();
    2996         else
    2997           value *= -CoinDrand48();
     3118            fabs(columnUpperWork_[iColumn])) {
     3119          value *= (1.0-m1 +m1*CoinDrand48());
     3120          value2 *= (1.0-m1 +m1*CoinDrand48());
     3121        } else {
     3122          value *= -(1.0-m1+m1*CoinDrand48());
     3123          value2 *= -(1.0-m1+m1*CoinDrand48());
     3124        }
    29983125      } else if (columnUpperWork_[iColumn]<largeValue_) {
    2999         value *= -CoinDrand48();
     3126        value *= -(1.0-m1+m1*CoinDrand48());
     3127        value2 *= -(1.0-m1+m1*CoinDrand48());
    30003128      } else {
    30013129        value=0.0;
    30023130      }
    3003       objectiveWork_[iColumn] += value;
    3004     }
    3005   }
     3131      if (value) {
     3132        int length = lengths[iColumn];
     3133        if (length>3) {
     3134          length = (int) ((double) length * factor);
     3135          length = max(3,length);
     3136        }
     3137        double multiplier;
     3138#if 1
     3139        if (length<10)
     3140          multiplier=weight[length];
     3141        else
     3142          multiplier = weight[10];
     3143#else
     3144        if (length<10)
     3145          multiplier=weight[length];
     3146        else
     3147          multiplier = weight[10]+extraWeight*(length-10);
     3148        multiplier *= 0.5;
     3149#endif
     3150        value *= multiplier;
     3151        value = min (value,value2);
     3152        if (fabs(value)<=dualTolerance_)
     3153          value=0.0;
     3154        if (currentValue) {
     3155          largest = max(largest,fabs(value));
     3156          if (fabs(value)>fabs(currentValue)*largestPerCent)
     3157            largestPerCent=fabs(value/currentValue);
     3158        } else {
     3159          largestZero=max(largestZero,fabs(value));
     3160        }
     3161        if (printOut)
     3162          printf("col %d cost %g change %g\n",iColumn,objectiveWork_[iColumn],value);
     3163        objectiveWork_[iColumn] += value;
     3164      }
     3165    }
     3166  }
     3167  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
     3168    <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
     3169    <<CoinMessageEol;
    30063170  // and zero changes
    30073171  //int nTotal = numberRows_+numberColumns_;
  • branches/pre/IdiSolve.cpp

    r63 r208  
    11021102  result.objval-=saveOffset;
    11031103  result.weighted=result.objval+weight*result.sumSquared;
    1104   if ((logLevel_&1)!=0) {
    1105     printf("%d Final infeas %g, obj %g - wtObj %g\n",
    1106                iter,result.infeas,result.objval,
    1107            result.weighted);
    1108   }
    11091104  return result;
    11101105}
  • branches/pre/Idiot.cpp

    r180 r208  
    1010#include "Idiot.hpp"
    1111#include "CoinHelperFunctions.hpp"
     12#include "CoinMessageHandler.hpp"
    1213// Redefine stuff for Clp
    1314#ifdef CLP_IDIOT
     15#include "ClpMessage.hpp"
    1416#define OsiObjOffset ClpObjOffset
    1517#endif
     
    105107}
    106108void
    107 Idiot::crash(int numberPass)
     109Idiot::crash(int numberPass, CoinMessageHandler * handler,const CoinMessages *messages)
    108110{
    109111  // lightweight options
     
    128130  else
    129131    majorIterations_=numberPass;
    130   printf("setting mu to %g and doing %d passes\n",mu_,majorIterations_);
    131   solve2();
     132  //printf("setting mu to %g and doing %d passes\n",mu_,majorIterations_);
     133  solve2(handler,messages);
    132134#ifdef CLP_IDIOT
    133135  crossOver(3); // make basis ?
     
    137139Idiot::solve()
    138140{
    139   solve2();
     141  CoinMessages dummy;
     142  solve2(NULL,&dummy);
    140143}
    141144void
    142 Idiot::solve2()
     145Idiot::solve2(CoinMessageHandler * handler,const CoinMessages * messages)
    143146{
    144147  int strategy=0;
    145148  double d2;
    146   double startTime = CoinCpuTime();
    147149  int i,n;
    148150  int allOnes=1;
     
    154156  int maxIts=maxIts_;
    155157  int logLevel=logLevel_;
     158  if (handler) {
     159    if (handler->logLevel()>0&&handler->logLevel()<3)
     160      logLevel=1;
     161    else if (!handler->logLevel())
     162      logLevel=0;
     163    else
     164      logLevel=7;
     165  }
    156166  double djExit=djTolerance_;
    157167  double djFlag=1.0+100.0*djExit;
     
    246256    ordEnd=ncols;
    247257  }
    248   if (offset) {
     258  if (offset&&logLevel>2) {
    249259    printf("** Objective offset is %g\n",offset);
    250260  }
     
    300310  if ((saveStrategy&8)!=0) strategy |= 64; /* don't allow large theta */
    301311  memcpy(saveSol,colsol,ncols*sizeof(double));
     312 
    302313  lastResult=IdiSolve(nrows,ncols,rowsol ,colsol,pi,
    303314                       dj,cost,rowlower,rowupper,
     
    305316                       0,mu,drop,
    306317                       maxmin,offset,strategy,djTol,djExit,djFlag);
     318  n=0;
     319  for (i=ordStart;i<ordEnd;i++) {
     320    if (colsol[i]>lower[i]+fixTolerance) {
     321      if (colsol[i]<upper[i]-fixTolerance) {
     322        n++;
     323      } else {
     324        colsol[i]=upper[i];
     325      }
     326      whenUsed[i]=iteration;
     327    } else {
     328      colsol[i]=lower[i];
     329    }
     330  }
     331  if ((logLevel_&1)!=0) {
     332#ifdef CLP_IDIOT
     333    if (!handler) {
     334#endif
     335      printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n",
     336             iteration,lastResult.infeas,lastResult.objval,mu,lastResult.iteration,n);
     337#ifdef CLP_IDIOT
     338    } else {
     339      handler->message(CLP_IDIOT_ITERATION,*messages)
     340        <<iteration<<lastResult.infeas<<lastResult.objval<<mu<<lastResult.iteration<<n
     341        <<CoinMessageEol;
     342    }
     343#endif
     344  }
    307345  iterationTotal = lastResult.iteration;
    308346  firstInfeas=lastResult.infeas;
     
    321359      bestFeasible=lastResult.objval;
    322360    }
     361    if ((saveStrategy&4096)) strategy |=256;
     362    if ((saveStrategy&4)!=0&&iteration>2) {
     363      /* go to relative tolerances */
     364      double weighted=10.0+fabs(lastWeighted);
     365      djExit=djTolerance_*weighted;
     366      djFlag=2.0*djExit;
     367      drop=drop_*weighted;
     368      djTol=0.01*djExit;
     369    }
     370    result=IdiSolve(nrows,ncols,rowsol ,colsol,pi,dj,
     371                     cost,rowlower,rowupper,
     372                     lower,upper,elemXX,row,columnStart,columnLength,lambda,
     373                     maxIts,mu,drop,
     374                     maxmin,offset,strategy,djTol,djExit,djFlag);
    323375    n=0;
    324376    for (i=ordStart;i<ordEnd;i++) {
     
    334386      }
    335387    }
    336     if ((logLevel&1)!=0) {
    337       printf(
    338               "Iteration %d - mu %g, drop %g, dj test %g, %d away from bounds\n",
    339               iteration,mu,drop,djExit,n);
    340     } else {
    341       printf(
    342               "%d - mu %g, infeasibility %g, objective %g, %d interior\n",
    343               iteration,mu,lastResult.infeas,lastResult.objval,n);
    344     }
    345     if ((saveStrategy&4096)) strategy |=256;
    346     if ((saveStrategy&4)!=0&&iteration>2) {
    347       /* go to relative tolerances */
    348       double weighted=10.0+fabs(lastWeighted);
    349       djExit=djTolerance_*weighted;
    350       djFlag=2.0*djExit;
    351       drop=drop_*weighted;
    352       djTol=0.01*djExit;
    353     }
    354     result=IdiSolve(nrows,ncols,rowsol ,colsol,pi,dj,
    355                      cost,rowlower,rowupper,
    356                      lower,upper,elemXX,row,columnStart,columnLength,lambda,
    357                      maxIts,mu,drop,
    358                      maxmin,offset,strategy,djTol,djExit,djFlag);
     388    if ((logLevel_&1)!=0) {
     389#ifdef CLP_IDIOT
     390      if (!handler) {
     391#endif
     392        printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n",
     393               iteration,result.infeas,result.objval,mu,result.iteration,n);
     394#ifdef CLP_IDIOT
     395      } else {
     396        handler->message(CLP_IDIOT_ITERATION,*messages)
     397          <<iteration<<result.infeas<<result.objval<<mu<<result.iteration<<n
     398          <<CoinMessageEol;
     399      }
     400#endif
     401    }
    359402    keepinfeas = result.infeas;
    360403    lastWeighted=result.weighted;
     
    492535    }
    493536    if (iteration>=majorIterations_) {
    494       printf("Exiting due to number of major iterations\n");
     537      if (logLevel>2)
     538        printf("Exiting due to number of major iterations\n");
    495539      break;
    496540    }
     
    560604      }
    561605    }
    562     printf("largest infeasibility is %g\n",large);
     606    if (logLevel>2)
     607      printf("largest infeasibility is %g\n",large);
    563608  }
    564609  /* subtract out lambda */
     
    603648  delete [] pi;
    604649  delete [] dj;
    605   printf("Idiot took %g seconds\n", CoinCpuTime()-startTime);
    606650  return ;
    607651}
  • branches/pre/Makefile.Clp

    r204 r208  
    11# Static or shared libraries should be built (STATIC or SHARED)?
    22LibType := SHARED
    3 #LibType := STATIC
     3LibType := STATIC
    44
    55# Select optimization (-O or -g). -O will be automatically bumped up to the
    66# highest level of optimization the compiler supports. If want something in
    77# between then specify the exact level you want, e.g., -O1 or -O2
    8 OptLevel := -g
    98OptLevel := -O3
    10 
     9#OptLevel := -g
     10#OptLevel := -O1
    1111
    1212LIBNAME := Clp
     
    3434LIBSRC += ClpSimplexPrimal.cpp
    3535LIBSRC += ClpSimplexPrimalQuadratic.cpp
    36 #LIBSRC += ClpPrimalQuadraticDantzig.cpp
     36LIBSRC += ClpSolve.cpp
    3737# and Presolve stuff
    3838LIBSRC += ClpPresolve.cpp
  • branches/pre/Test/ClpMain.cpp

    r207 r208  
    1414#include "CoinPragma.hpp"
    1515#include "CoinHelperFunctions.hpp"
    16 #define CLPVERSION "0.98.1"
     16#define CLPVERSION "0.98.5"
    1717
    1818#include "CoinMpsIO.hpp"
     
    2727#include "ClpPrimalColumnSteepest.hpp"
    2828#include "ClpPrimalColumnDantzig.hpp"
    29 
    3029#include "ClpPresolve.hpp"
    31 #ifdef CLP_IDIOT
    32 #include "Idiot.hpp"
    33 #endif
     30
    3431// For Branch and bound
    3532//  #include "OsiClpSolverInterface.hpp"
     
    5754  DUALTOLERANCE=1,PRIMALTOLERANCE,DUALBOUND,PRIMALWEIGHT,MAXTIME,OBJSCALE,
    5855
    59   LOGLEVEL=101,MAXFACTOR,PERTURBATION,MAXITERATION,PRESOLVEPASS,IDIOT,
     56  LOGLEVEL=101,MAXFACTOR,PERTVALUE,MAXITERATION,PRESOLVEPASS,IDIOT,
    6057 
    6158  DIRECTION=201,DUALPIVOT,SCALING,ERRORSALLOWED,KEEPNAMES,SPARSEFACTOR,
    62   PRIMALPIVOT,PRESOLVE,CRASH,BIASLU,
     59  PRIMALPIVOT,PRESOLVE,CRASH,BIASLU,PERTURBATION,
    6360 
    6461  DIRECTORY=301,IMPORT,EXPORT,RESTORE,SAVE,DUALSIMPLEX,PRIMALSIMPLEX,
     
    526523      model->setOptimizationDirection(value);
    527524      break;
    528     case PERTURBATION:
     525    case PERTVALUE:
    529526      model->setPerturbation(value);
    530527      break;
     
    560557    }
    561558    break;
    562   case PERTURBATION:
     559  case PERTVALUE:
    563560    value=model->perturbation();
    564561    break;
     
    575572#include <readline/readline.h>
    576573#include <readline/history.h>
    577 #include <signal.h>
    578 static ClpSimplex * currentModel = NULL;
    579 static void signal_handler(int whichSignal)
    580 {
    581   if (currentModel!=NULL)
    582     currentModel->setMaximumIterations(0); // stop at next iterations
    583   return;
    584 }
    585574#endif
    586575// Returns next valid field
     
    789778    parameters[numberParameters++]=
    790779      ClpItem("dualP!ivot","Dual pivot choice algorithm",
    791               "steep!est",DUALPIVOT);
     780              "auto!matic",DUALPIVOT);
    792781    parameters[numberParameters-1].append("dant!zig");
    793782    parameters[numberParameters-1].append("partial");
     783    parameters[numberParameters-1].append("steep!est");
    794784    parameters[numberParameters++]=
    795785      ClpItem("dualS!implex","Do dual simplex algorithm",
     
    857847              0,100,PRESOLVEPASS);
    858848    parameters[numberParameters++]=
    859       ClpItem("pert!urbation","Method of perturbation",
    860               -5000,102,PERTURBATION);
     849      ClpItem("pertV!alue","Method of perturbation",
     850              -5000,102,PERTVALUE,false);
     851    parameters[numberParameters++]=
     852      ClpItem("perturb!ation","Whether to perturb problem",
     853              "on",PERTURBATION);
     854    parameters[numberParameters-1].append("off");
    861855    parameters[numberParameters++]=
    862856      ClpItem("plus!Minus","Tries to make +- 1 matrix",
     
    866860              "on",PRESOLVE);
    867861    parameters[numberParameters-1].append("off");
     862    parameters[numberParameters-1].append("more");
    868863    parameters[numberParameters++]=
    869864      ClpItem("primalP!ivot","Primal pivot choice algorithm",
     
    929924    bool * goodModels = new bool[1];
    930925    int getNewMatrix=0;
    931 #ifdef READLINE     
    932     currentModel = models;
    933     // register signal handler
    934     signal(SIGINT,signal_handler);
    935 #endif
    936926   
    937927    // default action on import
     
    944934    // set reasonable defaults
    945935    int preSolve=5;
    946     models->setPerturbation(53);
     936    models->setPerturbation(50);
    947937    //models[0].scaling(1);
    948938    //models[0].setDualBound(1.0e6);
     
    955945    std::string directory ="./";
    956946    std::string field;
     947    std::cout<<"Coin LP version "<<CLPVERSION
     948             <<", build "<<__DATE__<<std::endl;
    957949   
    958950    while (1) {
     
    972964            <<"Enter ? for list of commands, (-)unitTest or (-)netlib"
    973965            <<" for tests"<<std::endl;
    974           break;
     966          field="-";
    975967        } else {
    976968          break;
     
    10841076            case DUALPIVOT:
    10851077              if (action==0) {
    1086                 ClpDualRowSteepest steep;
     1078                ClpDualRowSteepest steep(3);
    10871079                models[iModel].setDualRowPivotAlgorithm(steep);
    10881080              } else if (action==1) {
    10891081                ClpDualRowDantzig dantzig;
    10901082                models[iModel].setDualRowPivotAlgorithm(dantzig);
    1091               } else {
     1083              } else if (action==2) {
    10921084                // partial steep
    10931085                ClpDualRowSteepest steep(2);
     1086                models[iModel].setDualRowPivotAlgorithm(steep);
     1087              } else {
     1088                ClpDualRowSteepest steep;
    10941089                models[iModel].setDualRowPivotAlgorithm(steep);
    10951090              }
     
    11161111              models[iModel].factorization()->setBiasLU(action);
    11171112              break;
     1113            case PERTURBATION:
     1114              if (action==0)
     1115                models[iModel].setPerturbation(50);
     1116              else
     1117                models[iModel].setPerturbation(100);
     1118              break;
    11181119            case ERRORSALLOWED:
    11191120              allowImportErrors = action;
     
    11231124              break;
    11241125            case PRESOLVE:
    1125               preSolve = (1-action)*5;
     1126              if (action==0)
     1127                preSolve = 5;
     1128              else if (action==1)
     1129                preSolve=0;
     1130              else
     1131                preSolve=10;
    11261132              break;
    11271133            case CRASH:
     
    11401146          case PRIMALSIMPLEX:
    11411147            if (goodModels[iModel]) {
    1142               int saveMaxIterations = models[iModel].maximumIterations();
    1143               int finalStatus=-1;
    1144               int numberIterations=0;
    1145               time1 = CoinCpuTime();
    1146               ClpMatrixBase * saveMatrix=NULL;
     1148              ClpSimplex::SolveType method;
     1149              ClpSimplex::PresolveType presolveType;
    11471150              ClpSimplex * model2 = models+iModel;
    1148               ClpPresolve pinfo;
    1149               double timePresolve=0.0;
    1150               if (preSolve) {
    1151                 model2 = pinfo.presolvedModel(models[iModel],1.0e-8,
    1152                                               false,preSolve,true);
    1153                 timePresolve = CoinCpuTime()-time1;
    1154                
    1155                 std::cout<<"Presolve took "<<timePresolve<<" seconds"<<std::endl;
    1156                 if (model2) {
    1157                   //model2->checkSolution();
    1158                   if (type==DUALSIMPLEX) {
    1159                     int numberInfeasibilities = model2->tightenPrimalBounds();
    1160                     if (numberInfeasibilities) {
    1161                       std::cout<<"** Analysis indicates model infeasible"
    1162                                <<std::endl;
    1163                       model2 = models+iModel;
    1164                       preSolve=0;
    1165                     }
    1166                   }
    1167 #ifdef READLINE     
    1168                 currentModel = model2;
    1169 #endif
    1170                 } else {
    1171                   std::cout<<"** Analysis indicates model infeasible"
    1172                            <<std::endl;
    1173                   model2 = models+iModel;
    1174                   preSolve=0;
    1175                 }
    1176               }
    1177               if (getNewMatrix) {
    1178                 saveMatrix = model2->clpMatrix();
    1179                 ClpPackedMatrix* clpMatrix =
    1180                   dynamic_cast< ClpPackedMatrix*>(saveMatrix);
    1181                 if (clpMatrix) {
    1182                   if (getNewMatrix==1) {
    1183                     ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
    1184                     if (newMatrix->getIndices()) {
    1185                       std::cout<<"** Matrix is valid +- one"<<std::endl;
    1186                       model2->replaceMatrix(newMatrix);
    1187                     } else {
    1188                       std::cout<<"** Matrix is NOT valid +- one"<<std::endl;
    1189                       saveMatrix=NULL;
    1190                       delete newMatrix;
    1191                     }
    1192                   } else if (getNewMatrix==2) {
    1193                     ClpNetworkMatrix * newMatrix = new ClpNetworkMatrix(*(clpMatrix->matrix()));
    1194                     if (newMatrix->getIndices()) {
    1195                       std::cout<<"** Matrix is valid network"<<std::endl;
    1196                       model2->replaceMatrix(newMatrix);
    1197                     } else {
    1198                       std::cout<<"** Matrix is NOT valid network"<<std::endl;
    1199                       saveMatrix=NULL;
    1200                       delete newMatrix;
    1201                     }
    1202                   }
    1203                 } else {
    1204                   saveMatrix=NULL;
    1205                 }
    1206               }
    1207               if (model2->factorizationFrequency()==200) {
    1208                 // User did not touch preset
    1209                 model2->setFactorizationFrequency(100+model2->numberRows()/100);
    1210               }
    1211               if (type==DUALSIMPLEX) {
    1212                 if (doIdiot==-1)
    1213                   model2->crash(1000,1);
    1214                 //int status =model2->saveModel("xx.save");
    1215                 model2->dual();
    1216                 //int status =model2->saveModel("xx.save");
    1217               } else {
    1218 #ifdef CLP_IDIOT
    1219                 if (doIdiot==-2||doIdiot>0) {
    1220                   if (doIdiot==-2) {
    1221                     int numberColumns = model2->numberColumns();
    1222                     int numberRows = model2->numberRows();
    1223                     if (numberRows>2000&&numberColumns>2*numberRows) {
    1224                       ClpPackedMatrix* clpMatrix =
    1225                         dynamic_cast< ClpPackedMatrix*>(model2->clpMatrix());
    1226                       if (clpMatrix) {
    1227                         ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
    1228                         if (newMatrix->getIndices()) {
    1229                           saveMatrix = model2->clpMatrix();
    1230                           std::cout<<"** Matrix is valid +- one"<<std::endl;
    1231                           model2->replaceMatrix(newMatrix);
    1232                           int nPasses = 10+numberColumns/1000;
    1233                           nPasses = min(nPasses,100);
    1234                           Idiot info(*model2);
    1235                           info.crash(nPasses);
    1236                         } else {
    1237                           delete newMatrix;
    1238                           int nPasses = 10+numberColumns/100000;
    1239                           if (numberColumns>4*numberRows)
    1240                             nPasses = min(nPasses,50);
    1241                           else
    1242                             nPasses=5;
    1243                           Idiot info(*model2);
    1244                           info.crash(nPasses);
    1245                         }
    1246                       } else {
    1247                         ClpPlusMinusOneMatrix* clpMatrix =
    1248                           dynamic_cast< ClpPlusMinusOneMatrix*>(model2->clpMatrix());
    1249                         if (clpMatrix) {
    1250                           int nPasses = 10+numberColumns/1000;
    1251                           nPasses = min(nPasses,100);
    1252                           Idiot info(*model2);
    1253                           info.crash(nPasses);
    1254                         } else {
    1255                           int nPasses = 10+numberColumns/100000;
    1256                           if (numberColumns>4*numberRows)
    1257                             nPasses = min(nPasses,50);
    1258                           else
    1259                             nPasses=5;
    1260                           Idiot info(*model2);
    1261                           info.crash(nPasses);
    1262                         }
    1263                       }
    1264                     }
    1265                   } else {
    1266                     Idiot info(*model2);
    1267                     info.crash(doIdiot);
    1268                   }
    1269                 }
    1270 #endif
    1271                 int savePerturbation = model2->perturbation();
    1272                 if (savePerturbation==53)
    1273                   model2->setPerturbation(100);
    1274                 model2->primal(1);
    1275                 model2->setPerturbation(savePerturbation);
    1276               }
    1277               if (saveMatrix) {
    1278                 // delete and replace
    1279                 delete model2->clpMatrix();
    1280                 model2->replaceMatrix(saveMatrix);
    1281               }
    1282               numberIterations = model2->numberIterations();
    1283               finalStatus=model2->status();
    1284               if (preSolve) {
    1285                 double timeX=CoinCpuTime();
    1286                 pinfo.postsolve(true);
    1287                 timePresolve += CoinCpuTime()-timeX;
    1288 
    1289                 delete model2;
    1290 
    1291 #ifdef READLINE     
    1292                 currentModel = models+iModel;
    1293 #endif
    1294                 models[iModel].checkSolution();
    1295                 if (finalStatus&&finalStatus!=3) {
    1296                   printf("Resolving from postsolved model\n");
    1297                  
    1298                   int savePerturbation = models[iModel].perturbation();
    1299                   models[iModel].setPerturbation(100);
    1300                   models[iModel].primal(1);
    1301                   models[iModel].setPerturbation(savePerturbation);
    1302                   numberIterations += models[iModel].numberIterations();
    1303                   finalStatus=models[iModel].status();
    1304                 }
    1305               }
    1306               models[iModel].setMaximumIterations(saveMaxIterations);
    1307               time2 = CoinCpuTime();
    1308               totalTime += time2-time1;
    1309               std::cout<<"Result "<<finalStatus<<
    1310                 " - "<<models[iModel].objectiveValue()<<
    1311                 " iterations "<<numberIterations<<
    1312                 " took "<<time2-time1<<" seconds - total "<<totalTime;
    1313               if (preSolve)
    1314                 std::cout<<" (Presolve took "<<timePresolve<<")";
    1315               std::cout<<std::endl;
    1316               if (finalStatus)
    1317                 std::cerr<<"Non zero status "<<finalStatus<<
    1318                   std::endl;
    1319               time1=time2;
     1151              if (preSolve>5)
     1152                presolveType=ClpSimplex::presolveMaximum;
     1153              else if (preSolve)
     1154                presolveType=ClpSimplex::presolveOn;
     1155              else
     1156                presolveType=ClpSimplex::presolveOff;
     1157              if (type==DUALSIMPLEX)
     1158                method=ClpSimplex::useDual;
     1159              else
     1160                method=ClpSimplex::usePrimal;
     1161              int option=0;
     1162              if (model2->perturbation()==100||method==ClpSimplex::usePrimal)
     1163                option |= 1;
     1164              if (!model2->scalingFlag())
     1165                option |= 2;
     1166              if (doIdiot==-1)
     1167                option |=4;
     1168              if (doIdiot==0)
     1169                option |= 8;
     1170              model2->initialSolve(method,presolveType,option);
    13201171            } else {
    13211172              std::cout<<"** Current model not valid"<<std::endl;
     
    16611512            std::cout<<"Non default values:-"<<std::endl;
    16621513            std::cout<<"Perturbation "<<models[0].perturbation()
    1663                      <<" (default 100), Presolve being done with 5 passes"<<std::endl;
    1664             break;
     1514                     <<" (default 100), Presolve being done with 5 passes"
     1515                     <<std::endl;
     1516            std::cout <<"Dual steepest edge steep/partial on matrix shape"
     1517                      <<std::endl;
     1518            std::cout <<"If Factorization frequency default then done on size of matrix"
     1519                      <<std::endl;
     1520            break;
    16651521          case SOLUTION:
    16661522            if (goodModels[iModel]) {
  • branches/pre/Test/Makefile.test

    r204 r208  
    6363endif
    6464#LDFLAGS += -lefence
     65LDFLAGS += -static
     66CXXFLAGS += -pg
     67#if DENSE and using ATLAS
     68#LDFLAGS += -L/home/forrest/ATLAS/lib/Linux_P4SSE2 -llapack -lcblas -latlas -lf77blas -lg2c
     69#if DENSE and using given libraries
     70LDFLAGS += -llapack -lblas -lg2c
    6571
    6672###############################################################################
     
    8894        @mkdir -p $(TARGETDIR)
    8995        @rm -f $@
    90 #       $(CXX) -static $(CXXFLAGS) -o $@ $(TESTOBJ) $(LDFLAGS) $(SYSLD) -lm
    9196        $(CXX) $(CXXFLAGS) -o $@ $(TESTOBJ) $(LDFLAGS) $(SYSLD) -lm
    9297        ${CP} $@ ..
  • branches/pre/Test/unitTest.cpp

    r206 r208  
    288288          if (doIdiot>0) {
    289289            Idiot info(*model2);
    290             info.crash(doIdiot);
     290            info.crash(doIdiot,model2->messageHandler(),model2->messagesPointer());
    291291          }
    292292#endif
     
    347347          if (doIdiot>0) {
    348348            Idiot info(solution);
    349             info.crash(doIdiot);
     349            info.crash(doIdiot,solution.messageHandler(),solution.messagesPointer());
    350350          }
    351351#endif
  • branches/pre/include/ClpDualRowSteepest.hpp

    r205 r208  
    6262  //@{
    6363  /** Default Constructor
    64       0 is uninitialized, 1 full, 2 is partial uninitialized.
     64      0 is uninitialized, 1 full, 2 is partial uninitialized,
     65      3 starts as 2 but may switch to 1.
    6566      By partial is meant that the weights are updated as normal
    6667      but only part of the infeasible basic variables are scanned. 
    6768      This can be faster on very easy problems.
    6869  */
    69   ClpDualRowSteepest(int mode=0);
     70  ClpDualRowSteepest(int mode=3);
    7071 
    7172  /// Copy constructor
     
    100101  int state_;
    101102  /** If 0 then we are using uninitialized weights, 1 then full,
    102       if 2 then uninitialized partial */
     103      if 2 then uninitialized partial, 3 switchable */
    103104  int mode_;
    104105  /// weight array
  • branches/pre/include/ClpMessage.hpp

    r190 r208  
    6363  CLP_QUADRATIC_BOTH,
    6464  CLP_QUADRATIC_PRIMAL_DETAILS,
     65  CLP_IDIOT_ITERATION,
     66  CLP_INFEASIBLE,
     67  CLP_MATRIX_CHANGE,
     68  CLP_TIMING,
     69  CLP_INTERVAL_TIMING,
    6570  CLP_DUMMY_END
    6671};
  • branches/pre/include/ClpSimplex.hpp

    r205 r208  
    4949public:
    5050
     51  /** enums for solve function */
     52  enum SolveType {
     53    useDual=0,
     54    usePrimal,
     55    useSprint,
     56    automatic
     57  };
     58  enum PresolveType {
     59    presolveOn=0,
     60    presolveOff,
     61    presolveMaximum
     62  };
    5163  /** enums for status of various sorts.
    5264      First 4 match CoinWarmStartBasis,
     
    90102  /// Destructor
    91103   ~ClpSimplex (  );
    92   // Ones below are just ClpModel with setti
     104  // Ones below are just ClpModel with some changes
    93105  /** Loads a problem (the constraints on the
    94106        rows are given by lower and upper bounds). If a pointer is 0 then the
     
    143155  /**@name Functions most useful to user */
    144156  //@{
     157  /** General solve algorithm which can do presolve
     158      special options (bits)
     159      1 - do not perturb
     160      2 - do not scale
     161      4 - use crash (default allslack in dual, idiot in primal)
     162      8 - all slack basis in primal
     163      16 - switch off interrupt handling
     164      32 - do not try and make plus minus one matrix
     165   */
     166  int initialSolve(SolveType method=useDual, PresolveType presolve=presolveOn,
     167                   int specialOptions=0);
    145168  /** Dual algorithm - see ClpSimplexDual.hpp for method */
    146169  int dual(int ifValuesPass=0);
     
    157180  /// Passes in factorization
    158181  void setFactorization( ClpFactorization & factorization);
    159   /// Sets or unsets scaling, 0 -off, 1 on, 2 dynamic(later)
     182  /// Sets or unsets scaling, 0 -off, 1 equilibrium, 2 geometric, 3 dynamic(later)
    160183  void scaling(int mode=1);
    161184  /// Gets scalingFlag
     
    856879  /// Column scale factors
    857880  double * columnScale_;
    858   /// Scale flag, 0 none, 1 normal, 2 dynamic
     881  /// Scale flag, 0 none, 1 equilibrium, 2 geometric, 3 dynamic
    859882  int scalingFlag_;
    860883  /// Number of times code has tentatively thought optimal
  • branches/pre/include/Idiot.hpp

    r61 r208  
    1010typedef int CoinBigIndex;
    1111#endif
     12class CoinMessageHandler;
     13class CoinMessages;
    1214// for use internally
    1315typedef struct {
     
    6567  void solve();
    6668  /// Lightweight "crash"
    67   void crash(int numberPass=0);
     69  void crash(int numberPass,CoinMessageHandler * handler ,const CoinMessages * messages);
    6870  /** Use simplex to get an optimal solution
    6971      mode is how many steps the simplex crossover should take to
     
    149151
    150152  /// Does actual work
    151   void solve2();
     153  void solve2(CoinMessageHandler * handler,const CoinMessages *messages);
    152154IdiotResult IdiSolve(
    153155                     int nrows, int ncols, double * rowsol , double * colsol,
Note: See TracChangeset for help on using the changeset viewer.