266 | | // Set the lb/ub on the duals |
267 | | volprob.dsize = dsize; |
268 | | volprob.psize = psize; |
269 | | volprob.dual_lb.allocate(dsize); |
270 | | volprob.dual_ub.allocate(dsize); |
271 | | int i; |
272 | | const double * rowLower = model->rowLower(); |
273 | | const double * rowUpper = model->rowUpper(); |
274 | | for (i = 0; i < dsize; ++i) { |
275 | | double range; |
276 | | convertBoundToSense(rowLower[i], rowUpper[i], |
277 | | sense[i], rhs[i], range); |
278 | | switch (sense[i]) { |
279 | | case 'E': |
280 | | volprob.dual_lb[i] = -1.0e31; |
281 | | volprob.dual_ub[i] = 1.0e31; |
282 | | break; |
283 | | case 'L': |
284 | | volprob.dual_lb[i] = -1.0e31; |
285 | | volprob.dual_ub[i] = 0.0; |
286 | | break; |
287 | | case 'G': |
288 | | volprob.dual_lb[i] = 0.0; |
289 | | volprob.dual_ub[i] = 1.0e31; |
290 | | break; |
291 | | default: |
292 | | printf("Volume Algorithm can't work if there is a non ELG row\n"); |
293 | | return 1; |
294 | | } |
295 | | } |
296 | | // Check bounds |
297 | | double * saveLower = model->columnLower(); |
298 | | double * saveUpper = model->columnUpper(); |
299 | | bool good = true; |
300 | | for (i = 0; i < psize; i++) { |
301 | | if (saveLower[i] < -1.0e20 || saveUpper[i] > 1.0e20) { |
302 | | good = false; |
303 | | break; |
304 | | } |
305 | | } |
306 | | if (!good) { |
307 | | saveLower = CoinCopyOfArray(model->columnLower(), psize); |
308 | | saveUpper = CoinCopyOfArray(model->columnUpper(), psize); |
309 | | for (i = 0; i < psize; i++) { |
310 | | if (saveLower[i] < -1.0e20) |
311 | | saveLower[i] = -1.0e20; |
312 | | if(saveUpper[i] > 1.0e20) |
313 | | saveUpper[i] = 1.0e20; |
314 | | } |
315 | | } |
316 | | lpHook myHook(saveLower, saveUpper, |
317 | | model->objective(), |
318 | | rhs, sense, *mat); |
| 266 | // Set the lb/ub on the duals |
| 267 | volprob.dsize = dsize; |
| 268 | volprob.psize = psize; |
| 269 | volprob.dual_lb.allocate(dsize); |
| 270 | volprob.dual_ub.allocate(dsize); |
| 271 | int i; |
| 272 | const double *rowLower = model->rowLower(); |
| 273 | const double *rowUpper = model->rowUpper(); |
| 274 | for (i = 0; i < dsize; ++i) { |
| 275 | double range; |
| 276 | convertBoundToSense(rowLower[i], rowUpper[i], |
| 277 | sense[i], rhs[i], range); |
| 278 | switch (sense[i]) { |
| 279 | case 'E': |
| 280 | volprob.dual_lb[i] = -1.0e31; |
| 281 | volprob.dual_ub[i] = 1.0e31; |
| 282 | break; |
| 283 | case 'L': |
| 284 | volprob.dual_lb[i] = -1.0e31; |
| 285 | volprob.dual_ub[i] = 0.0; |
| 286 | break; |
| 287 | case 'G': |
| 288 | volprob.dual_lb[i] = 0.0; |
| 289 | volprob.dual_ub[i] = 1.0e31; |
| 290 | break; |
| 291 | default: |
| 292 | printf("Volume Algorithm can't work if there is a non ELG row\n"); |
| 293 | return 1; |
| 294 | } |
| 295 | } |
| 296 | // Check bounds |
| 297 | double *saveLower = model->columnLower(); |
| 298 | double *saveUpper = model->columnUpper(); |
| 299 | bool good = true; |
| 300 | for (i = 0; i < psize; i++) { |
| 301 | if (saveLower[i] < -1.0e20 || saveUpper[i] > 1.0e20) { |
| 302 | good = false; |
| 303 | break; |
| 304 | } |
| 305 | } |
| 306 | if (!good) { |
| 307 | saveLower = CoinCopyOfArray(model->columnLower(), psize); |
| 308 | saveUpper = CoinCopyOfArray(model->columnUpper(), psize); |
| 309 | for (i = 0; i < psize; i++) { |
| 310 | if (saveLower[i] < -1.0e20) |
| 311 | saveLower[i] = -1.0e20; |
| 312 | if (saveUpper[i] > 1.0e20) |
| 313 | saveUpper[i] = 1.0e20; |
| 314 | } |
| 315 | } |
| 316 | lpHook myHook(saveLower, saveUpper, |
| 317 | model->objective(), |
| 318 | rhs, sense, *mat); |
865 | | debugInt[0]=numberRows(); |
866 | | debugInt[1]=numberColumns(); |
867 | | debugInt[2]=matrix()->getNumElements(); |
868 | | #endif |
869 | | if (!objective_ || !matrix_) { |
870 | | // totally empty |
871 | | handler_->message(CLP_EMPTY_PROBLEM, messages_) |
872 | | << 0 |
873 | | << 0 |
874 | | << 0 |
875 | | << CoinMessageEol; |
876 | | return -1; |
877 | | } else if (!numberRows_ || !numberColumns_ || !getNumElements()) { |
878 | | presolve = ClpSolve::presolveOff; |
879 | | } |
880 | | if (objective_->type() >= 2 && optimizationDirection_ == 0) { |
881 | | // pretend linear |
882 | | savedObjective = objective_; |
883 | | // make up objective |
884 | | double * obj = new double[numberColumns_]; |
885 | | for (int i = 0; i < numberColumns_; i++) { |
886 | | double l = fabs(columnLower_[i]); |
887 | | double u = fabs(columnUpper_[i]); |
888 | | obj[i] = 0.0; |
889 | | if (CoinMin(l, u) < 1.0e20) { |
890 | | if (l < u) |
891 | | obj[i] = 1.0 + randomNumberGenerator_.randomDouble() * 1.0e-2; |
892 | | else |
893 | | obj[i] = -1.0 - randomNumberGenerator_.randomDouble() * 1.0e-2; |
894 | | } |
895 | | } |
896 | | objective_ = new ClpLinearObjective(obj, numberColumns_); |
897 | | delete [] obj; |
898 | | } |
899 | | ClpSimplex * model2 = this; |
900 | | bool interrupt = (options.getSpecialOption(2) == 0); |
901 | | CoinSighandler_t saveSignal = static_cast<CoinSighandler_t> (0); |
902 | | if (interrupt) { |
903 | | currentModel = model2; |
904 | | // register signal handler |
905 | | saveSignal = signal(SIGINT, signal_handler); |
906 | | } |
907 | | // If no status array - set up basis |
908 | | if (!status_) |
909 | | allSlackBasis(); |
910 | | ClpPresolve * pinfo = new ClpPresolve(); |
911 | | pinfo->setSubstitution(options.substitution()); |
912 | | int presolveOptions = options.presolveActions(); |
913 | | bool presolveToFile = (presolveOptions & 0x40000000) != 0; |
914 | | presolveOptions &= ~0x40000000; |
915 | | if ((presolveOptions & 0xffffff) != 0) |
916 | | pinfo->setPresolveActions(presolveOptions); |
917 | | // switch off singletons to slacks |
918 | | //pinfo->setDoSingletonColumn(false); // done by bits |
919 | | int printOptions = options.getSpecialOption(5); |
920 | | if ((printOptions & 1) != 0) |
921 | | pinfo->statistics(); |
922 | | double timePresolve = 0.0; |
923 | | double timeIdiot = 0.0; |
924 | | double timeCore = 0.0; |
925 | | eventHandler()->event(ClpEventHandler::presolveStart); |
926 | | int savePerturbation = perturbation_; |
927 | | int saveScaling = scalingFlag_; |
| 865 | debugInt[0] = numberRows(); |
| 866 | debugInt[1] = numberColumns(); |
| 867 | debugInt[2] = matrix()->getNumElements(); |
| 868 | #endif |
| 869 | if (!objective_ || !matrix_) { |
| 870 | // totally empty |
| 871 | handler_->message(CLP_EMPTY_PROBLEM, messages_) |
| 872 | << 0 |
| 873 | << 0 |
| 874 | << 0 |
| 875 | << CoinMessageEol; |
| 876 | return -1; |
| 877 | } else if (!numberRows_ || !numberColumns_ || !getNumElements()) { |
| 878 | presolve = ClpSolve::presolveOff; |
| 879 | } |
| 880 | if (objective_->type() >= 2 && optimizationDirection_ == 0) { |
| 881 | // pretend linear |
| 882 | savedObjective = objective_; |
| 883 | // make up objective |
| 884 | double *obj = new double[numberColumns_]; |
| 885 | for (int i = 0; i < numberColumns_; i++) { |
| 886 | double l = fabs(columnLower_[i]); |
| 887 | double u = fabs(columnUpper_[i]); |
| 888 | obj[i] = 0.0; |
| 889 | if (CoinMin(l, u) < 1.0e20) { |
| 890 | if (l < u) |
| 891 | obj[i] = 1.0 + randomNumberGenerator_.randomDouble() * 1.0e-2; |
| 892 | else |
| 893 | obj[i] = -1.0 - randomNumberGenerator_.randomDouble() * 1.0e-2; |
| 894 | } |
| 895 | } |
| 896 | objective_ = new ClpLinearObjective(obj, numberColumns_); |
| 897 | delete[] obj; |
| 898 | } |
| 899 | ClpSimplex *model2 = this; |
| 900 | bool interrupt = (options.getSpecialOption(2) == 0); |
| 901 | CoinSighandler_t saveSignal = static_cast< CoinSighandler_t >(0); |
| 902 | if (interrupt) { |
| 903 | currentModel = model2; |
| 904 | // register signal handler |
| 905 | saveSignal = signal(SIGINT, signal_handler); |
| 906 | } |
| 907 | // If no status array - set up basis |
| 908 | if (!status_) |
| 909 | allSlackBasis(); |
| 910 | ClpPresolve *pinfo = new ClpPresolve(); |
| 911 | pinfo->setSubstitution(options.substitution()); |
| 912 | int presolveOptions = options.presolveActions(); |
| 913 | bool presolveToFile = (presolveOptions & 0x40000000) != 0; |
| 914 | presolveOptions &= ~0x40000000; |
| 915 | if ((presolveOptions & 0xffffff) != 0) |
| 916 | pinfo->setPresolveActions(presolveOptions); |
| 917 | // switch off singletons to slacks |
| 918 | //pinfo->setDoSingletonColumn(false); // done by bits |
| 919 | int printOptions = options.getSpecialOption(5); |
| 920 | if ((printOptions & 1) != 0) |
| 921 | pinfo->statistics(); |
| 922 | double timePresolve = 0.0; |
| 923 | double timeIdiot = 0.0; |
| 924 | double timeCore = 0.0; |
| 925 | eventHandler()->event(ClpEventHandler::presolveStart); |
| 926 | int savePerturbation = perturbation_; |
| 927 | int saveScaling = scalingFlag_; |
1016 | | if (model2->primalTolerance()==1.0e-7&&model2->dualTolerance()==1.0e-7) { |
1017 | | model2->setPrimalTolerance(CLP_NEW_TOLERANCE); |
1018 | | model2->setDualTolerance(CLP_NEW_TOLERANCE); |
1019 | | } |
1020 | | } |
1021 | | #endif |
1022 | | #endif |
1023 | | model2->eventHandler()->setSimplex(model2); |
1024 | | int rcode=model2->eventHandler()->event(ClpEventHandler::presolveSize); |
1025 | | // see if too big or small |
1026 | | if (rcode==2) { |
1027 | | delete model2; |
1028 | | delete pinfo; |
1029 | | return -2; |
1030 | | } else if (rcode==3) { |
1031 | | delete model2; |
1032 | | delete pinfo; |
1033 | | return -3; |
1034 | | } |
| 1016 | if (model2->primalTolerance() == 1.0e-7 && model2->dualTolerance() == 1.0e-7) { |
| 1017 | model2->setPrimalTolerance(CLP_NEW_TOLERANCE); |
| 1018 | model2->setDualTolerance(CLP_NEW_TOLERANCE); |
| 1019 | } |
| 1020 | } |
| 1021 | #endif |
| 1022 | #endif |
| 1023 | model2->eventHandler()->setSimplex(model2); |
| 1024 | int rcode = model2->eventHandler()->event(ClpEventHandler::presolveSize); |
| 1025 | // see if too big or small |
| 1026 | if (rcode == 2) { |
| 1027 | delete model2; |
| 1028 | delete pinfo; |
| 1029 | return -2; |
| 1030 | } else if (rcode == 3) { |
| 1031 | delete model2; |
| 1032 | delete pinfo; |
| 1033 | return -3; |
| 1034 | } |
| 1035 | } |
| 1036 | model2->setMoreSpecialOptions(model2->moreSpecialOptions() & (~1024)); |
| 1037 | model2->eventHandler()->setSimplex(model2); |
| 1038 | // We may be better off using original (but if dual leave because of bounds) |
| 1039 | if (presolve != ClpSolve::presolveOff && numberRows_ < 1.01 * model2->numberRows_ && numberColumns_ < 1.01 * model2->numberColumns_ |
| 1040 | && model2 != this) { |
| 1041 | if (method != ClpSolve::useDual || (numberRows_ == model2->numberRows_ && numberColumns_ == model2->numberColumns_)) { |
| 1042 | delete model2; |
| 1043 | model2 = this; |
| 1044 | presolve = ClpSolve::presolveOff; |
| 1045 | } |
| 1046 | } |
| 1047 | } |
| 1048 | #ifdef CLP_USEFUL_PRINTOUT |
| 1049 | debugInt[3] = model2->numberRows(); |
| 1050 | debugInt[4] = model2->numberColumns(); |
| 1051 | debugInt[5] = model2->matrix()->getNumElements(); |
| 1052 | // analyze |
| 1053 | { |
| 1054 | double time1 = CoinCpuTime(); |
| 1055 | int numberColumns = model2->numberColumns(); |
| 1056 | const double *columnLower = model2->columnLower(); |
| 1057 | const double *columnUpper = model2->columnUpper(); |
| 1058 | int numberRows = model2->numberRows(); |
| 1059 | const double *rowLower = model2->rowLower(); |
| 1060 | const double *rowUpper = model2->rowUpper(); |
| 1061 | const double *objective = model2->objective(); |
| 1062 | CoinPackedMatrix *matrix = model2->matrix(); |
| 1063 | CoinBigIndex numberElements = matrix->getNumElements(); |
| 1064 | const int *columnLength = matrix->getVectorLengths(); |
| 1065 | //const CoinBigIndex * columnStart = matrix->getVectorStarts(); |
| 1066 | const double *elementByColumn = matrix->getElements(); |
| 1067 | const int *row = matrix->getIndices(); |
| 1068 | int *rowCount = new int[numberRows]; |
| 1069 | memset(rowCount, 0, numberRows * sizeof(int)); |
| 1070 | int n = CoinMax(2 * numberRows, numberElements); |
| 1071 | n = CoinMax(2 * numberColumns, n); |
| 1072 | double *check = new double[n]; |
| 1073 | memcpy(check, elementByColumn, numberElements * sizeof(double)); |
| 1074 | for (int i = 0; i < numberElements; i++) { |
| 1075 | check[i] = fabs(check[i]); |
| 1076 | rowCount[row[i]]++; |
| 1077 | } |
| 1078 | int largestIndex = 0; |
| 1079 | for (int i = 0; i < numberColumns; i++) { |
| 1080 | largestIndex = CoinMax(largestIndex, columnLength[i]); |
| 1081 | } |
| 1082 | debugInt[12] = largestIndex; |
| 1083 | largestIndex = 0; |
| 1084 | for (int i = 0; i < numberRows; i++) { |
| 1085 | largestIndex = CoinMax(largestIndex, rowCount[i]); |
| 1086 | } |
| 1087 | n = numberElements; |
| 1088 | delete[] rowCount; |
| 1089 | debugInt[11] = largestIndex; |
| 1090 | std::sort(check, check + n); |
| 1091 | debugDouble[4] = check[0]; |
| 1092 | debugDouble[5] = check[n - 1]; |
| 1093 | int nAtOne = 0; |
| 1094 | int nDifferent = 0; |
| 1095 | double last = -COIN_DBL_MAX; |
| 1096 | for (int i = 0; i < n; i++) { |
| 1097 | if (fabs(last - check[i]) > 1.0e-12) { |
| 1098 | nDifferent++; |
| 1099 | last = check[i]; |
| 1100 | } |
| 1101 | if (check[i] == 1.0) |
| 1102 | nAtOne++; |
| 1103 | } |
| 1104 | debugInt[10] = nDifferent; |
| 1105 | debugInt[15] = nAtOne; |
| 1106 | int nInf = 0; |
| 1107 | int nZero = 0; |
| 1108 | n = 0; |
| 1109 | for (int i = 0; i < numberRows; i++) { |
| 1110 | double value = fabs(rowLower[i]); |
| 1111 | if (!value) |
| 1112 | nZero++; |
| 1113 | else if (value != COIN_DBL_MAX) |
| 1114 | check[n++] = value; |
| 1115 | else |
| 1116 | nInf++; |
| 1117 | } |
| 1118 | for (int i = 0; i < numberRows; i++) { |
| 1119 | double value = fabs(rowUpper[i]); |
| 1120 | if (!value) |
| 1121 | nZero++; |
| 1122 | else if (value != COIN_DBL_MAX) |
| 1123 | check[n++] = value; |
| 1124 | else |
| 1125 | nInf++; |
| 1126 | } |
| 1127 | debugInt[16] = nInf; |
| 1128 | debugInt[20] = nZero; |
| 1129 | if (n) { |
| 1130 | std::sort(check, check + n); |
| 1131 | debugDouble[0] = check[0]; |
| 1132 | debugDouble[1] = check[n - 1]; |
| 1133 | } else { |
| 1134 | debugDouble[0] = 0.0; |
| 1135 | debugDouble[1] = 0.0; |
| 1136 | } |
| 1137 | nAtOne = 0; |
| 1138 | nDifferent = 0; |
| 1139 | last = -COIN_DBL_MAX; |
| 1140 | for (int i = 0; i < n; i++) { |
| 1141 | if (fabs(last - check[i]) > 1.0e-12) { |
| 1142 | nDifferent++; |
| 1143 | last = check[i]; |
| 1144 | } |
| 1145 | if (check[i] == 1.0) |
| 1146 | nAtOne++; |
| 1147 | } |
| 1148 | debugInt[8] = nDifferent; |
| 1149 | debugInt[13] = nAtOne; |
| 1150 | nZero = 0; |
| 1151 | n = 0; |
| 1152 | for (int i = 0; i < numberColumns; i++) { |
| 1153 | double value = fabs(objective[i]); |
| 1154 | if (value) |
| 1155 | check[n++] = value; |
| 1156 | else |
| 1157 | nZero++; |
| 1158 | } |
| 1159 | debugInt[21] = nZero; |
| 1160 | if (n) { |
| 1161 | std::sort(check, check + n); |
| 1162 | debugDouble[2] = check[0]; |
| 1163 | debugDouble[3] = check[n - 1]; |
| 1164 | } else { |
| 1165 | debugDouble[2] = 0.0; |
| 1166 | debugDouble[3] = 0.0; |
| 1167 | } |
| 1168 | nAtOne = 0; |
| 1169 | nDifferent = 0; |
| 1170 | last = -COIN_DBL_MAX; |
| 1171 | for (int i = 0; i < n; i++) { |
| 1172 | if (fabs(last - check[i]) > 1.0e-12) { |
| 1173 | nDifferent++; |
| 1174 | last = check[i]; |
| 1175 | } |
| 1176 | if (check[i] == 1.0) |
| 1177 | nAtOne++; |
| 1178 | } |
| 1179 | debugInt[9] = nDifferent; |
| 1180 | debugInt[14] = nAtOne; |
| 1181 | nInf = 0; |
| 1182 | nZero = 0; |
| 1183 | n = 0; |
| 1184 | for (int i = 0; i < numberColumns; i++) { |
| 1185 | double value = fabs(columnLower[i]); |
| 1186 | if (!value) |
| 1187 | nZero++; |
| 1188 | else if (value != COIN_DBL_MAX) |
| 1189 | check[n++] = value; |
| 1190 | else |
| 1191 | nInf++; |
| 1192 | } |
| 1193 | for (int i = 0; i < numberColumns; i++) { |
| 1194 | double value = fabs(columnUpper[i]); |
| 1195 | if (!value) |
| 1196 | nZero++; |
| 1197 | else if (value != COIN_DBL_MAX) |
| 1198 | check[n++] = value; |
| 1199 | else |
| 1200 | nInf++; |
| 1201 | } |
| 1202 | debugInt[17] = nInf; |
| 1203 | double smallestColBound; |
| 1204 | double largestColBound; |
| 1205 | if (n) { |
| 1206 | std::sort(check, check + n); |
| 1207 | smallestColBound = check[0]; |
| 1208 | largestColBound = check[n - 1]; |
| 1209 | } else { |
| 1210 | smallestColBound = 0.0; |
| 1211 | largestColBound = 0.0; |
| 1212 | } |
| 1213 | nAtOne = 0; |
| 1214 | nDifferent = 0; |
| 1215 | last = -COIN_DBL_MAX; |
| 1216 | for (int i = 0; i < n; i++) { |
| 1217 | if (fabs(last - check[i]) > 1.0e-12) { |
| 1218 | nDifferent++; |
| 1219 | last = check[i]; |
| 1220 | } |
| 1221 | if (check[i] == 1.0) |
| 1222 | nAtOne++; |
| 1223 | } |
| 1224 | //debugInt[8]=nDifferent; |
| 1225 | //debugInt[13]=nAtOne; |
| 1226 | printf("BENCHMARK_STATS rhs %d different - %g -> %g (%d at one, %d infinite, %d zero)\n", |
| 1227 | debugInt[8], debugDouble[0], debugDouble[1], debugInt[13], debugInt[16], debugInt[20]); |
| 1228 | printf("BENCHMARK_STATS col %d different - %g -> %g (%d at one, %d infinite, %d zero)\n", |
| 1229 | nDifferent, smallestColBound, largestColBound, nAtOne, nInf, nZero); |
| 1230 | printf("BENCHMARK_STATS els %d different - %g -> %g (%d at one) - longest r,c %d,%d\n", |
| 1231 | debugInt[10], debugDouble[4], debugDouble[5], debugInt[15], |
| 1232 | debugInt[11], debugInt[12]); |
| 1233 | printf("BENCHMARK_STATS obj %d different - %g -> %g (%d at one, %d zero) - time %g\n", |
| 1234 | debugInt[9], debugDouble[2], debugDouble[3], debugInt[14], debugInt[21], |
| 1235 | CoinCpuTime() - time1); |
| 1236 | delete[] check; |
| 1237 | } |
| 1238 | #endif |
| 1239 | if (interrupt) |
| 1240 | currentModel = model2; |
| 1241 | int saveMoreOptions = moreSpecialOptions_; |
| 1242 | // For below >0 overrides |
| 1243 | // 0 means no, -1 means maybe |
| 1244 | int doIdiot = 0; |
| 1245 | int doCrash = 0; |
| 1246 | int doSprint = 0; |
| 1247 | int doSlp = 0; |
| 1248 | int primalStartup = 1; |
| 1249 | model2->eventHandler()->event(ClpEventHandler::presolveBeforeSolve); |
| 1250 | #if CLP_POOL_MATRIX |
| 1251 | if (vectorMode() >= 10) { |
| 1252 | ClpPoolMatrix *poolMatrix = new ClpPoolMatrix(*model2->matrix()); |
| 1253 | char output[80]; |
| 1254 | int numberDifferent = poolMatrix->getNumDifferentElements(); |
| 1255 | if (numberDifferent > 0) { |
| 1256 | sprintf(output, "Pool matrix has %d different values", |
| 1257 | numberDifferent); |
| 1258 | model2->replaceMatrix(poolMatrix, true); |
| 1259 | } else { |
| 1260 | delete poolMatrix; |
| 1261 | sprintf(output, "Pool matrix has more than %d different values - no good", |
| 1262 | -numberDifferent); |
| 1263 | } |
| 1264 | handler_->message(CLP_GENERAL, messages_) << output |
| 1265 | << CoinMessageEol; |
| 1266 | } |
| 1267 | #endif |
| 1268 | int tryItSave = 0; |
| 1269 | #if CLPSOLVE_ACTIONS |
| 1270 | if (method == ClpSolve::automatic) |
| 1271 | model2->moreSpecialOptions_ |= 8192; // stop switch over |
| 1272 | #endif |
| 1273 | // switch to primal from automatic if just one cost entry |
| 1274 | if (method == ClpSolve::automatic && model2->numberColumns() > 5000 |
| 1275 | #ifndef CLPSOLVE_ACTIONS |
| 1276 | && (specialOptions_ & 1024) != 0 |
| 1277 | #endif |
| 1278 | ) { |
| 1279 | // look at original model for objective |
| 1280 | int numberColumns = model2->numberColumns(); |
| 1281 | int numberRows = model2->numberRows(); |
| 1282 | int numberColumnsOrig = this->numberColumns(); |
| 1283 | const double *obj = this->objective(); |
| 1284 | int nNon = 0; |
| 1285 | bool allOnes = true; |
| 1286 | for (int i = 0; i < numberColumnsOrig; i++) { |
| 1287 | if (obj[i]) { |
| 1288 | nNon++; |
| 1289 | if (fabs(obj[i]) != 1.0) |
| 1290 | allOnes = false; |
| 1291 | } |
| 1292 | } |
| 1293 | if (nNon <= 1 || allOnes || (options.getExtraInfo(1) > 0 && options.getSpecialOption(1) == 2)) { |
| 1294 | #ifdef COIN_DEVELOP |
| 1295 | printf("Forcing primal\n"); |
| 1296 | #endif |
| 1297 | method = ClpSolve::usePrimal; |
| 1298 | #ifndef CLPSOLVE_ACTIONS |
| 1299 | tryItSave = (numberRows > 200 && numberColumns > 2000 && (numberColumns > 2 * numberRows || (specialOptions_ & 1024) != 0)) ? 3 : 0; |
| 1300 | #else |
| 1301 | if (numberRows > 200 && numberColumns > 2000 && numberColumns > 2 * numberRows) { |
| 1302 | tryItSave = 3; |
| 1303 | } else { |
| 1304 | // If rhs also rubbish then maybe |
| 1305 | int numberRowsOrig = this->numberRows(); |
| 1306 | const double *rowLower = this->rowLower(); |
| 1307 | const double *rowUpper = this->rowUpper(); |
| 1308 | double last = COIN_DBL_MAX; |
| 1309 | int nDifferent = 0; |
| 1310 | for (int i = 0; i < numberRowsOrig; i++) { |
| 1311 | double value = fabs(rowLower[i]); |
| 1312 | if (value && value != COIN_DBL_MAX) { |
| 1313 | if (value != last) { |
| 1314 | nDifferent++; |
| 1315 | last = value; |
| 1316 | } |
1049 | | } |
1050 | | #ifdef CLP_USEFUL_PRINTOUT |
1051 | | debugInt[3]=model2->numberRows(); |
1052 | | debugInt[4]=model2->numberColumns(); |
1053 | | debugInt[5]=model2->matrix()->getNumElements(); |
1054 | | // analyze |
1055 | | { |
1056 | | double time1 =CoinCpuTime(); |
1057 | | int numberColumns = model2->numberColumns(); |
1058 | | const double * columnLower = model2->columnLower(); |
1059 | | const double * columnUpper = model2->columnUpper(); |
1060 | | int numberRows = model2->numberRows(); |
1061 | | const double * rowLower = model2->rowLower(); |
1062 | | const double * rowUpper = model2->rowUpper(); |
1063 | | const double * objective = model2->objective(); |
1064 | | CoinPackedMatrix * matrix = model2->matrix(); |
1065 | | CoinBigIndex numberElements = matrix->getNumElements(); |
1066 | | const int * columnLength = matrix->getVectorLengths(); |
1067 | | //const CoinBigIndex * columnStart = matrix->getVectorStarts(); |
1068 | | const double * elementByColumn = matrix->getElements(); |
1069 | | const int * row = matrix->getIndices(); |
1070 | | int * rowCount = new int [numberRows]; |
1071 | | memset(rowCount,0,numberRows*sizeof(int)); |
1072 | | int n=CoinMax (2*numberRows,numberElements); |
1073 | | n= CoinMax(2*numberColumns,n); |
1074 | | double * check = new double[n]; |
1075 | | memcpy(check,elementByColumn,numberElements*sizeof(double)); |
1076 | | for (int i=0;i<numberElements;i++) { |
1077 | | check[i]=fabs(check[i]); |
1078 | | rowCount[row[i]]++; |
1079 | | } |
1080 | | int largestIndex=0; |
1081 | | for (int i=0;i<numberColumns;i++) { |
1082 | | largestIndex=CoinMax(largestIndex,columnLength[i]); |
1083 | | } |
1084 | | debugInt[12]=largestIndex; |
1085 | | largestIndex=0; |
1086 | | for (int i=0;i<numberRows;i++) { |
1087 | | largestIndex=CoinMax(largestIndex,rowCount[i]); |
1088 | | } |
1089 | | n=numberElements; |
1090 | | delete [] rowCount; |
1091 | | debugInt[11]=largestIndex; |
1092 | | std::sort(check,check+n); |
1093 | | debugDouble[4]=check[0]; |
1094 | | debugDouble[5]=check[n-1]; |
1095 | | int nAtOne=0; |
1096 | | int nDifferent=0; |
1097 | | double last=-COIN_DBL_MAX; |
1098 | | for (int i=0;i<n;i++) { |
1099 | | if (fabs(last-check[i])>1.0e-12) { |
1100 | | nDifferent++; |
1101 | | last=check[i]; |
1102 | | } |
1103 | | if (check[i]==1.0) |
1104 | | nAtOne++; |
1105 | | } |
1106 | | debugInt[10]=nDifferent; |
1107 | | debugInt[15]=nAtOne; |
1108 | | int nInf=0; |
1109 | | int nZero=0; |
1110 | | n=0; |
1111 | | for (int i=0;i<numberRows;i++) { |
1112 | | double value=fabs(rowLower[i]); |
1113 | | if (!value) |
1114 | | nZero++; |
1115 | | else if (value!=COIN_DBL_MAX) |
1116 | | check[n++]=value; |
1117 | | else |
1118 | | nInf++; |
1119 | | } |
1120 | | for (int i=0;i<numberRows;i++) { |
1121 | | double value=fabs(rowUpper[i]); |
1122 | | if (!value) |
1123 | | nZero++; |
1124 | | else if (value!=COIN_DBL_MAX) |
1125 | | check[n++]=value; |
1126 | | else |
1127 | | nInf++; |
1128 | | } |
1129 | | debugInt[16]=nInf; |
1130 | | debugInt[20]=nZero; |
1131 | | if (n) { |
1132 | | std::sort(check,check+n); |
1133 | | debugDouble[0]=check[0]; |
1134 | | debugDouble[1]=check[n-1]; |
1135 | | } else { |
1136 | | debugDouble[0]=0.0; |
1137 | | debugDouble[1]=0.0; |
1138 | | } |
1139 | | nAtOne=0; |
1140 | | nDifferent=0; |
1141 | | last=-COIN_DBL_MAX; |
1142 | | for (int i=0;i<n;i++) { |
1143 | | if (fabs(last-check[i])>1.0e-12) { |
1144 | | nDifferent++; |
1145 | | last=check[i]; |
1146 | | } |
1147 | | if (check[i]==1.0) |
1148 | | nAtOne++; |
1149 | | } |
1150 | | debugInt[8]=nDifferent; |
1151 | | debugInt[13]=nAtOne; |
1152 | | nZero=0; |
1153 | | n=0; |
1154 | | for (int i=0;i<numberColumns;i++) { |
1155 | | double value=fabs(objective[i]); |
1156 | | if (value) |
1157 | | check[n++]=value; |
1158 | | else |
1159 | | nZero++; |
1160 | | } |
1161 | | debugInt[21]=nZero; |
1162 | | if (n) { |
1163 | | std::sort(check,check+n); |
1164 | | debugDouble[2]=check[0]; |
1165 | | debugDouble[3]=check[n-1]; |
1166 | | } else { |
1167 | | debugDouble[2]=0.0; |
1168 | | debugDouble[3]=0.0; |
1169 | | } |
1170 | | nAtOne=0; |
1171 | | nDifferent=0; |
1172 | | last=-COIN_DBL_MAX; |
1173 | | for (int i=0;i<n;i++) { |
1174 | | if (fabs(last-check[i])>1.0e-12) { |
1175 | | nDifferent++; |
1176 | | last=check[i]; |
1177 | | } |
1178 | | if (check[i]==1.0) |
1179 | | nAtOne++; |
1180 | | } |
1181 | | debugInt[9]=nDifferent; |
1182 | | debugInt[14]=nAtOne; |
1183 | | nInf=0; |
1184 | | nZero=0; |
1185 | | n=0; |
1186 | | for (int i=0;i<numberColumns;i++) { |
1187 | | double value=fabs(columnLower[i]); |
1188 | | if (!value) |
1189 | | nZero++; |
1190 | | else if (value!=COIN_DBL_MAX) |
1191 | | check[n++]=value; |
1192 | | else |
1193 | | nInf++; |
1194 | | } |
1195 | | for (int i=0;i<numberColumns;i++) { |
1196 | | double value=fabs(columnUpper[i]); |
1197 | | if (!value) |
1198 | | nZero++; |
1199 | | else if (value!=COIN_DBL_MAX) |
1200 | | check[n++]=value; |
1201 | | else |
1202 | | nInf++; |
1203 | | } |
1204 | | debugInt[17]=nInf; |
1205 | | double smallestColBound; |
1206 | | double largestColBound; |
1207 | | if (n) { |
1208 | | std::sort(check,check+n); |
1209 | | smallestColBound=check[0]; |
1210 | | largestColBound=check[n-1]; |
1211 | | } else { |
1212 | | smallestColBound=0.0; |
1213 | | largestColBound=0.0; |
1214 | | } |
1215 | | nAtOne=0; |
1216 | | nDifferent=0; |
1217 | | last=-COIN_DBL_MAX; |
1218 | | for (int i=0;i<n;i++) { |
1219 | | if (fabs(last-check[i])>1.0e-12) { |
1220 | | nDifferent++; |
1221 | | last=check[i]; |
1222 | | } |
1223 | | if (check[i]==1.0) |
1224 | | nAtOne++; |
1225 | | } |
1226 | | //debugInt[8]=nDifferent; |
1227 | | //debugInt[13]=nAtOne; |
1228 | | printf("BENCHMARK_STATS rhs %d different - %g -> %g (%d at one, %d infinite, %d zero)\n", |
1229 | | debugInt[8],debugDouble[0],debugDouble[1],debugInt[13],debugInt[16],debugInt[20]); |
1230 | | printf("BENCHMARK_STATS col %d different - %g -> %g (%d at one, %d infinite, %d zero)\n", |
1231 | | nDifferent,smallestColBound,largestColBound,nAtOne,nInf,nZero); |
1232 | | printf("BENCHMARK_STATS els %d different - %g -> %g (%d at one) - longest r,c %d,%d\n", |
1233 | | debugInt[10],debugDouble[4],debugDouble[5],debugInt[15], |
1234 | | debugInt[11],debugInt[12]); |
1235 | | printf("BENCHMARK_STATS obj %d different - %g -> %g (%d at one, %d zero) - time %g\n", |
1236 | | debugInt[9],debugDouble[2],debugDouble[3],debugInt[14],debugInt[21], |
1237 | | CoinCpuTime()-time1); |
1238 | | delete [] check; |
1239 | | } |
1240 | | #endif |
1241 | | if (interrupt) |
1242 | | currentModel = model2; |
1243 | | int saveMoreOptions = moreSpecialOptions_; |
1244 | | // For below >0 overrides |
1245 | | // 0 means no, -1 means maybe |
1246 | | int doIdiot = 0; |
1247 | | int doCrash = 0; |
1248 | | int doSprint = 0; |
1249 | | int doSlp = 0; |
1250 | | int primalStartup = 1; |
1251 | | model2->eventHandler()->event(ClpEventHandler::presolveBeforeSolve); |
1252 | | #if CLP_POOL_MATRIX |
1253 | | if (vectorMode()>=10) { |
1254 | | ClpPoolMatrix * poolMatrix = new ClpPoolMatrix(*model2->matrix()); |
1255 | | char output[80]; |
1256 | | int numberDifferent=poolMatrix->getNumDifferentElements(); |
1257 | | if (numberDifferent>0) { |
1258 | | sprintf(output,"Pool matrix has %d different values", |
1259 | | numberDifferent); |
1260 | | model2->replaceMatrix(poolMatrix,true); |
1261 | | } else { |
1262 | | delete poolMatrix; |
1263 | | sprintf(output,"Pool matrix has more than %d different values - no good", |
1264 | | -numberDifferent); |
1265 | | } |
1266 | | handler_->message(CLP_GENERAL, messages_) << output |
1267 | | << CoinMessageEol; |
1268 | | } |
1269 | | #endif |
1270 | | int tryItSave = 0; |
1271 | | #if CLPSOLVE_ACTIONS |
1272 | | if (method == ClpSolve::automatic ) |
1273 | | model2->moreSpecialOptions_ |= 8192; // stop switch over |
1274 | | #endif |
1275 | | // switch to primal from automatic if just one cost entry |
1276 | | if (method == ClpSolve::automatic && model2->numberColumns() > 5000 |
1277 | | #ifndef CLPSOLVE_ACTIONS |
1278 | | && (specialOptions_ & 1024) != 0 |
1279 | | #endif |
1280 | | ) { |
1281 | | // look at original model for objective |
1282 | | int numberColumns = model2->numberColumns(); |
1283 | | int numberRows = model2->numberRows(); |
1284 | | int numberColumnsOrig = this->numberColumns(); |
1285 | | const double * obj = this->objective(); |
1286 | | int nNon = 0; |
1287 | | bool allOnes=true; |
1288 | | for (int i = 0; i < numberColumnsOrig; i++) { |
1289 | | if (obj[i]) { |
1290 | | nNon++; |
1291 | | if (fabs(obj[i])!=1.0) |
1292 | | allOnes=false; |
1293 | | } |
1294 | | } |
1295 | | if (nNon <= 1 || allOnes || |
1296 | | (options.getExtraInfo(1) > 0 && options.getSpecialOption(1)==2)) { |
1297 | | #ifdef COIN_DEVELOP |
1298 | | printf("Forcing primal\n"); |
1299 | | #endif |
1300 | | method = ClpSolve::usePrimal; |
1301 | | #ifndef CLPSOLVE_ACTIONS |
1302 | | tryItSave = (numberRows > 200 && numberColumns > 2000 && |
1303 | | (numberColumns > 2 * numberRows || (specialOptions_ & 1024) != 0)) ? 3 : 0; |
| 1325 | } |
| 1326 | if (nDifferent < 2) |
| 1327 | tryItSave = 1; |
| 1328 | } |
| 1329 | #endif |
| 1330 | } |
| 1331 | } |
| 1332 | if (method != ClpSolve::useDual && method != ClpSolve::useBarrier |
| 1333 | && method != ClpSolve::tryBenders && method != ClpSolve::tryDantzigWolfe |
| 1334 | && method != ClpSolve::useBarrierNoCross) { |
| 1335 | switch (options.getSpecialOption(1)) { |
| 1336 | case 0: |
| 1337 | doIdiot = -1; |
| 1338 | doCrash = -1; |
| 1339 | doSprint = -1; |
| 1340 | break; |
| 1341 | case 1: |
| 1342 | doIdiot = 0; |
| 1343 | doCrash = 1; |
| 1344 | if (options.getExtraInfo(1) > 0) |
| 1345 | doCrash = options.getExtraInfo(1); |
| 1346 | doSprint = 0; |
| 1347 | break; |
| 1348 | case 2: |
| 1349 | doIdiot = 1; |
| 1350 | if (options.getExtraInfo(1) > 0) |
| 1351 | doIdiot = options.getExtraInfo(1); |
| 1352 | doCrash = 0; |
| 1353 | doSprint = 0; |
| 1354 | break; |
| 1355 | case 3: |
| 1356 | doIdiot = 0; |
| 1357 | doCrash = 0; |
| 1358 | doSprint = 1; |
| 1359 | break; |
| 1360 | case 4: |
| 1361 | doIdiot = 0; |
| 1362 | doCrash = 0; |
| 1363 | doSprint = 0; |
| 1364 | break; |
| 1365 | case 5: |
| 1366 | doIdiot = 0; |
| 1367 | doCrash = -1; |
| 1368 | doSprint = -1; |
| 1369 | break; |
| 1370 | case 6: |
| 1371 | doIdiot = -1; |
| 1372 | doCrash = -1; |
| 1373 | doSprint = 0; |
| 1374 | break; |
| 1375 | case 7: |
| 1376 | doIdiot = -1; |
| 1377 | doCrash = 0; |
| 1378 | doSprint = -1; |
| 1379 | break; |
| 1380 | case 8: |
| 1381 | doIdiot = -1; |
| 1382 | doCrash = 0; |
| 1383 | doSprint = 0; |
| 1384 | break; |
| 1385 | case 9: |
| 1386 | doIdiot = 0; |
| 1387 | doCrash = 0; |
| 1388 | doSprint = -1; |
| 1389 | break; |
| 1390 | case 10: |
| 1391 | doIdiot = 0; |
| 1392 | doCrash = 0; |
| 1393 | doSprint = 0; |
| 1394 | if (options.getExtraInfo(1)) |
| 1395 | doSlp = options.getExtraInfo(1); |
| 1396 | break; |
| 1397 | case 11: |
| 1398 | doIdiot = 0; |
| 1399 | doCrash = 0; |
| 1400 | doSprint = 0; |
| 1401 | primalStartup = 0; |
| 1402 | break; |
| 1403 | default: |
| 1404 | abort(); |
| 1405 | } |
| 1406 | } else if (method != ClpSolve::tryBenders && method != ClpSolve::tryDantzigWolfe) { |
| 1407 | // Dual |
| 1408 | switch (options.getSpecialOption(0)) { |
| 1409 | case 0: |
| 1410 | doIdiot = 0; |
| 1411 | doCrash = 0; |
| 1412 | doSprint = 0; |
| 1413 | break; |
| 1414 | case 1: |
| 1415 | doIdiot = 0; |
| 1416 | doCrash = 1; |
| 1417 | if (options.getExtraInfo(0) > 0) |
| 1418 | doCrash = options.getExtraInfo(0); |
| 1419 | doSprint = 0; |
| 1420 | break; |
| 1421 | case 2: |
| 1422 | doIdiot = -1; |
| 1423 | if (options.getExtraInfo(0) > 0) |
| 1424 | doIdiot = options.getExtraInfo(0); |
| 1425 | doCrash = 0; |
| 1426 | doSprint = 0; |
| 1427 | break; |
| 1428 | default: |
| 1429 | abort(); |
| 1430 | } |
| 1431 | } else { |
| 1432 | // decomposition |
| 1433 | } |
| 1434 | #ifndef NO_RTTI |
| 1435 | ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(objectiveAsObject())); |
1305 | | if (numberRows > 200 && numberColumns > 2000 && |
1306 | | numberColumns > 2 * numberRows) { |
1307 | | tryItSave= 3; |
1308 | | } else { |
1309 | | // If rhs also rubbish then maybe |
1310 | | int numberRowsOrig = this->numberRows(); |
1311 | | const double * rowLower = this->rowLower(); |
1312 | | const double * rowUpper = this->rowUpper(); |
1313 | | double last=COIN_DBL_MAX; |
1314 | | int nDifferent=0; |
1315 | | for (int i=0;i<numberRowsOrig;i++) { |
1316 | | double value=fabs(rowLower[i]); |
1317 | | if (value&&value!=COIN_DBL_MAX) { |
1318 | | if (value!=last) { |
1319 | | nDifferent++; |
1320 | | last=value; |
1321 | | } |
1322 | | } |
1323 | | value=fabs(rowUpper[i]); |
1324 | | if (value&&value!=COIN_DBL_MAX) { |
1325 | | if (value!=last) { |
1326 | | nDifferent++; |
1327 | | last=value; |
1328 | | } |
1329 | | } |
1330 | | } |
1331 | | if (nDifferent<2) |
1332 | | tryItSave=1; |
1333 | | } |
1334 | | #endif |
1335 | | } |
1336 | | } |
1337 | | if (method != ClpSolve::useDual && method != ClpSolve::useBarrier |
1338 | | && method != ClpSolve::tryBenders && method != ClpSolve::tryDantzigWolfe |
1339 | | && method != ClpSolve::useBarrierNoCross) { |
1340 | | switch (options.getSpecialOption(1)) { |
1341 | | case 0: |
1342 | | doIdiot = -1; |
1343 | | doCrash = -1; |
1344 | | doSprint = -1; |
1345 | | break; |
1346 | | case 1: |
1347 | | doIdiot = 0; |
1348 | | doCrash = 1; |
1349 | | if (options.getExtraInfo(1) > 0) |
1350 | | doCrash = options.getExtraInfo(1); |
1351 | | doSprint = 0; |
1352 | | break; |
1353 | | case 2: |
1354 | | doIdiot = 1; |
1355 | | if (options.getExtraInfo(1) > 0) |
1356 | | doIdiot = options.getExtraInfo(1); |
1357 | | doCrash = 0; |
1358 | | doSprint = 0; |
1359 | | break; |
1360 | | case 3: |
1361 | | doIdiot = 0; |
1362 | | doCrash = 0; |
1363 | | doSprint = 1; |
1364 | | break; |
1365 | | case 4: |
1366 | | doIdiot = 0; |
1367 | | doCrash = 0; |
1368 | | doSprint = 0; |
1369 | | break; |
1370 | | case 5: |
1371 | | doIdiot = 0; |
1372 | | doCrash = -1; |
1373 | | doSprint = -1; |
1374 | | break; |
1375 | | case 6: |
1376 | | doIdiot = -1; |
1377 | | doCrash = -1; |
1378 | | doSprint = 0; |
1379 | | break; |
1380 | | case 7: |
1381 | | doIdiot = -1; |
1382 | | doCrash = 0; |
1383 | | doSprint = -1; |
1384 | | break; |
1385 | | case 8: |
1386 | | doIdiot = -1; |
1387 | | doCrash = 0; |
1388 | | doSprint = 0; |
1389 | | break; |
1390 | | case 9: |
1391 | | doIdiot = 0; |
1392 | | doCrash = 0; |
1393 | | doSprint = -1; |
1394 | | break; |
1395 | | case 10: |
1396 | | doIdiot = 0; |
1397 | | doCrash = 0; |
1398 | | doSprint = 0; |
1399 | | if (options.getExtraInfo(1) ) |
1400 | | doSlp = options.getExtraInfo(1); |
1401 | | break; |
1402 | | case 11: |
1403 | | doIdiot = 0; |
1404 | | doCrash = 0; |
1405 | | doSprint = 0; |
1406 | | primalStartup = 0; |
1407 | | break; |
1408 | | default: |
1409 | | abort(); |
1410 | | } |
1411 | | } else if (method != ClpSolve::tryBenders && method != ClpSolve::tryDantzigWolfe) { |
1412 | | // Dual |
1413 | | switch (options.getSpecialOption(0)) { |
1414 | | case 0: |
1415 | | doIdiot = 0; |
1416 | | doCrash = 0; |
1417 | | doSprint = 0; |
1418 | | break; |
1419 | | case 1: |
1420 | | doIdiot = 0; |
1421 | | doCrash = 1; |
1422 | | if (options.getExtraInfo(0) > 0) |
1423 | | doCrash = options.getExtraInfo(0); |
1424 | | doSprint = 0; |
1425 | | break; |
1426 | | case 2: |
1427 | | doIdiot = -1; |
1428 | | if (options.getExtraInfo(0) > 0) |
1429 | | doIdiot = options.getExtraInfo(0); |
1430 | | doCrash = 0; |
1431 | | doSprint = 0; |
1432 | | break; |
1433 | | default: |
1434 | | abort(); |
1435 | | } |
1436 | | } else { |
1437 | | // decomposition |
1438 | | } |
1439 | | #ifndef NO_RTTI |
1440 | | ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objectiveAsObject())); |
1441 | | #else |
1442 | | ClpQuadraticObjective * quadraticObj = NULL; |
1443 | | if (objective_->type() == 2) |
1444 | | quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_)); |
1445 | | #endif |
1446 | | // If quadratic then primal or barrier or slp |
1447 | | if (quadraticObj) { |
1448 | | doSprint = 0; |
1449 | | //doIdiot = 0; |
1450 | | // off |
1451 | | if (method == ClpSolve::useBarrier) |
1452 | | method = ClpSolve::useBarrierNoCross; |
1453 | | else if (method != ClpSolve::useBarrierNoCross) |
1454 | | method = ClpSolve::usePrimal; |
1455 | | } |
| 1437 | ClpQuadraticObjective *quadraticObj = NULL; |
| 1438 | if (objective_->type() == 2) |
| 1439 | quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_)); |
| 1440 | #endif |
| 1441 | // If quadratic then primal or barrier or slp |
| 1442 | if (quadraticObj) { |
| 1443 | doSprint = 0; |
| 1444 | //doIdiot = 0; |
| 1445 | // off |
| 1446 | if (method == ClpSolve::useBarrier) |
| 1447 | method = ClpSolve::useBarrierNoCross; |
| 1448 | else if (method != ClpSolve::useBarrierNoCross) |
| 1449 | method = ClpSolve::usePrimal; |
| 1450 | } |
1467 | | if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) { |
1468 | | // network - switch off stuff |
| 1462 | if (dynamic_cast< ClpNetworkMatrix * >(matrix_)) { |
| 1463 | // network - switch off stuff |
| 1464 | doIdiot = 0; |
| 1465 | if (doSprint < 0) |
| 1466 | doSprint = 0; |
| 1467 | } |
| 1468 | #else |
| 1469 | if (matrix_->type() == 11) { |
| 1470 | // network - switch off stuff |
| 1471 | doIdiot = 0; |
| 1472 | //doSprint=0; |
| 1473 | } |
| 1474 | #endif |
| 1475 | #endif |
| 1476 | int numberColumns = model2->numberColumns(); |
| 1477 | int numberRows = model2->numberRows(); |
| 1478 | // If not all slack basis - switch off all except sprint |
| 1479 | int numberRowsBasic = 0; |
| 1480 | int iRow; |
| 1481 | for (iRow = 0; iRow < numberRows; iRow++) |
| 1482 | if (model2->getRowStatus(iRow) == basic) |
| 1483 | numberRowsBasic++; |
| 1484 | if (numberRowsBasic < numberRows && objective_->type() < 2) { |
| 1485 | doIdiot = 0; |
| 1486 | doCrash = 0; |
| 1487 | //doSprint=0; |
| 1488 | } |
| 1489 | if (options.getSpecialOption(3) == 0) { |
| 1490 | if (numberElements > 100000) |
| 1491 | plusMinus = true; |
| 1492 | if (numberElements > 10000 && (doIdiot || doSprint)) |
| 1493 | plusMinus = true; |
| 1494 | } else if ((specialOptions_ & 1024) != 0) { |
| 1495 | plusMinus = true; |
| 1496 | } |
| 1497 | #ifndef SLIM_CLP |
| 1498 | // Statistics (+1,-1, other) - used to decide on strategy if not +-1 |
| 1499 | CoinBigIndex statistics[3] = { -1, 0, 0 }; |
| 1500 | if (plusMinus) { |
| 1501 | saveMatrix = model2->clpMatrix(); |
| 1502 | #ifndef NO_RTTI |
| 1503 | ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(saveMatrix); |
| 1504 | #else |
| 1505 | ClpPackedMatrix *clpMatrix = NULL; |
| 1506 | if (saveMatrix->type() == 1) |
| 1507 | clpMatrix = static_cast< ClpPackedMatrix * >(saveMatrix); |
| 1508 | #endif |
| 1509 | if (clpMatrix) { |
| 1510 | ClpPlusMinusOneMatrix *newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix())); |
| 1511 | if (newMatrix->getIndices()) { |
| 1512 | // CHECKME This test of specialOptions and the one above |
| 1513 | // don't seem compatible. |
| 1514 | #ifndef ABC_INHERIT |
| 1515 | if ((specialOptions_ & 1024) == 0) { |
| 1516 | model2->replaceMatrix(newMatrix); |
| 1517 | } else { |
| 1518 | #endif |
| 1519 | // in integer (or abc) - just use for sprint/idiot |
| 1520 | saveMatrix = NULL; |
| 1521 | delete newMatrix; |
| 1522 | #ifndef ABC_INHERIT |
| 1523 | } |
| 1524 | #endif |
| 1525 | } else { |
| 1526 | handler_->message(CLP_MATRIX_CHANGE, messages_) |
| 1527 | << "+- 1" |
| 1528 | << CoinMessageEol; |
| 1529 | CoinMemcpyN(newMatrix->startPositive(), 3, statistics); |
| 1530 | saveMatrix = NULL; |
| 1531 | plusMinus = false; |
| 1532 | delete newMatrix; |
| 1533 | } |
| 1534 | } else { |
| 1535 | saveMatrix = NULL; |
| 1536 | plusMinus = false; |
| 1537 | } |
| 1538 | } |
| 1539 | #endif |
| 1540 | if (this->factorizationFrequency() == 200) { |
| 1541 | // User did not touch preset |
| 1542 | model2->defaultFactorizationFrequency(); |
| 1543 | } else if (model2 != this) { |
| 1544 | // make sure model2 has correct value |
| 1545 | model2->setFactorizationFrequency(this->factorizationFrequency()); |
| 1546 | } |
| 1547 | if (method == ClpSolve::automatic) { |
| 1548 | if (doSprint == 0 && doIdiot == 0) { |
| 1549 | // off |
| 1550 | method = ClpSolve::useDual; |
| 1551 | } else { |
| 1552 | // only do primal if sprint or idiot |
| 1553 | if (doSprint > 0) { |
| 1554 | method = ClpSolve::usePrimalorSprint; |
| 1555 | } else if (doIdiot > 0) { |
| 1556 | method = ClpSolve::usePrimal; |
| 1557 | } else { |
| 1558 | if (numberElements < 500000) { |
| 1559 | // Small problem |
| 1560 | if (numberRows * 10 > numberColumns || numberColumns < 6000 |
| 1561 | || (numberRows * 20 > numberColumns && !plusMinus)) |
| 1562 | doSprint = 0; // switch off sprint |
| 1563 | } else { |
| 1564 | // larger problem |
| 1565 | if (numberRows * 8 > numberColumns) |
| 1566 | doSprint = 0; // switch off sprint |
| 1567 | } |
| 1568 | // switch off idiot or sprint if any free variable |
| 1569 | // switch off sprint if very few with costs |
| 1570 | // or great variation in cost |
| 1571 | int iColumn; |
| 1572 | const double *columnLower = model2->columnLower(); |
| 1573 | const double *columnUpper = model2->columnUpper(); |
| 1574 | const double *objective = model2->objective(); |
| 1575 | int nObj = 0; |
| 1576 | int nFree = 0; |
| 1577 | double smallestObj = COIN_DBL_MAX; |
| 1578 | double largestObj = 0.0; |
| 1579 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
| 1580 | if (columnLower[iColumn] < -1.0e10 && columnUpper[iColumn] > 1.0e10) { |
| 1581 | nFree++; |
| 1582 | } else if (objective[iColumn]) { |
| 1583 | nObj++; |
| 1584 | smallestObj = CoinMin(smallestObj, objective[iColumn]); |
| 1585 | largestObj = CoinMax(largestObj, objective[iColumn]); |
| 1586 | } |
| 1587 | } |
| 1588 | if (nObj * 10 < numberColumns || smallestObj * 10.0 < largestObj) |
| 1589 | doSprint = 0; |
| 1590 | if (nFree) |
1470 | | if (doSprint < 0) |
1471 | | doSprint = 0; |
1472 | | } |
| 1592 | if (nFree * 10 > numberRows) |
| 1593 | doSprint = 0; |
| 1594 | int nPasses = 0; |
| 1595 | // look at rhs |
| 1596 | int iRow; |
| 1597 | double largest = 0.0; |
| 1598 | double smallest = 1.0e30; |
| 1599 | double largestGap = 0.0; |
| 1600 | int numberNotE = 0; |
| 1601 | bool notInteger = false; |
| 1602 | for (iRow = 0; iRow < numberRows; iRow++) { |
| 1603 | double value1 = model2->rowLower_[iRow]; |
| 1604 | if (value1 && value1 > -1.0e31) { |
| 1605 | largest = CoinMax(largest, fabs(value1)); |
| 1606 | smallest = CoinMin(smallest, fabs(value1)); |
| 1607 | if (fabs(value1 - floor(value1 + 0.5)) > 1.0e-8) { |
| 1608 | notInteger = true; |
| 1609 | break; |
| 1610 | } |
| 1611 | } |
| 1612 | double value2 = model2->rowUpper_[iRow]; |
| 1613 | if (value2 && value2 < 1.0e31) { |
| 1614 | largest = CoinMax(largest, fabs(value2)); |
| 1615 | smallest = CoinMin(smallest, fabs(value2)); |
| 1616 | if (fabs(value2 - floor(value2 + 0.5)) > 1.0e-8) { |
| 1617 | notInteger = true; |
| 1618 | break; |
| 1619 | } |
| 1620 | } |
| 1621 | // CHECKME This next bit can't be right... |
| 1622 | if (value2 > value1) { |
| 1623 | numberNotE++; |
| 1624 | //if (value2 > 1.0e31 || value1 < -1.0e31) |
| 1625 | // largestGap = COIN_DBL_MAX; |
| 1626 | //else |
| 1627 | // largestGap = value2 - value1; |
| 1628 | } |
| 1629 | } |
| 1630 | int tryIt = (numberRows > 200 && numberColumns > 2000 && (numberColumns > 2 * numberRows || (method != ClpSolve::useDual && (specialOptions_ & 1024) != 0))) ? 3 : 0; |
| 1631 | tryItSave = tryIt; |
| 1632 | if (numberRows < 1000 && numberColumns < 3000) |
| 1633 | tryIt = 0; |
| 1634 | if (notInteger) |
| 1635 | tryIt = 0; |
| 1636 | if (largest / smallest > 10 || (largest / smallest > 2.0 && largest > 50)) |
| 1637 | tryIt = 0; |
| 1638 | if (tryIt) { |
| 1639 | if (largest / smallest > 2.0) { |
| 1640 | nPasses = 10 + numberColumns / 100000; |
| 1641 | nPasses = CoinMin(nPasses, 50); |
| 1642 | nPasses = CoinMax(nPasses, 15); |
| 1643 | if (numberRows > 20000 && nPasses > 5) { |
| 1644 | // Might as well go for it |
| 1645 | nPasses = CoinMax(nPasses, 71); |
| 1646 | } else if (numberRows > 2000 && nPasses > 5) { |
| 1647 | nPasses = CoinMax(nPasses, 50); |
| 1648 | } else if (numberElements < 3 * numberColumns) { |
| 1649 | nPasses = CoinMin(nPasses, 10); // probably not worh it |
| 1650 | } |
| 1651 | } else if (largest / smallest > 1.01 || numberElements <= 3 * numberColumns) { |
| 1652 | nPasses = 10 + numberColumns / 1000; |
| 1653 | nPasses = CoinMin(nPasses, 100); |
| 1654 | nPasses = CoinMax(nPasses, 30); |
| 1655 | if (numberRows > 25000) { |
| 1656 | // Might as well go for it |
| 1657 | nPasses = CoinMax(nPasses, 71); |
| 1658 | } |
| 1659 | if (!largestGap) |
| 1660 | nPasses *= 2; |
| 1661 | } else { |
| 1662 | nPasses = 10 + numberColumns / 1000; |
| 1663 | nPasses = CoinMax(nPasses, 100); |
| 1664 | if (!largestGap) |
| 1665 | nPasses *= 2; |
| 1666 | nPasses = CoinMin(nPasses, 200); |
| 1667 | } |
| 1668 | } |
| 1669 | //printf("%d rows %d cols plus %c tryIt %c largest %g smallest %g largestGap %g npasses %d sprint %c\n", |
| 1670 | // numberRows,numberColumns,plusMinus ? 'Y' : 'N', |
| 1671 | // tryIt ? 'Y' :'N',largest,smallest,largestGap,nPasses,doSprint ? 'Y' :'N'); |
| 1672 | //exit(0); |
| 1673 | if (!tryIt || nPasses <= 5) |
| 1674 | doIdiot = 0; |
| 1675 | #if CLPSOLVE_ACTIONS |
| 1676 | if (doIdiot && doSprint < 0 && wasAutomatic && 20 * model2->numberRows() > model2->numberColumns()) |
| 1677 | doSprint = 0; // switch off sprint |
| 1678 | #endif |
| 1679 | if (doSprint) { |
| 1680 | method = ClpSolve::usePrimalorSprint; |
| 1681 | } else if (doIdiot) { |
| 1682 | method = ClpSolve::usePrimal; |
| 1683 | } else { |
| 1684 | method = ClpSolve::useDual; |
| 1685 | } |
| 1686 | } |
| 1687 | } |
| 1688 | } |
| 1689 | if (method == ClpSolve::tryBenders) { |
| 1690 | // Now build model |
| 1691 | int lengthNames = model2->lengthNames(); |
| 1692 | model2->setLengthNames(0); |
| 1693 | CoinModel *build = model2->createCoinModel(); |
| 1694 | model2->setLengthNames(lengthNames); |
| 1695 | CoinStructuredModel benders; |
| 1696 | build->convertMatrix(); |
| 1697 | int numberBlocks = options.independentOption(0); |
| 1698 | benders.setMessageHandler(handler_); |
| 1699 | numberBlocks = benders.decompose(*build, 2, numberBlocks, NULL); |
| 1700 | delete build; |
| 1701 | //exit(0); |
| 1702 | if (numberBlocks) { |
| 1703 | options.setIndependentOption(1, 1); // don't do final clean up |
| 1704 | model2->solveBenders(&benders, options); |
| 1705 | //move solution |
| 1706 | method = ClpSolve::notImplemented; |
| 1707 | time2 = CoinCpuTime(); |
| 1708 | timeCore = time2 - timeX; |
| 1709 | handler_->message(CLP_INTERVAL_TIMING, messages_) |
| 1710 | << "Crossover" << timeCore << time2 - time1 |
| 1711 | << CoinMessageEol; |
| 1712 | timeX = time2; |
| 1713 | } else { |
| 1714 | printf("No structure\n"); |
| 1715 | method = ClpSolve::useDual; |
| 1716 | } |
| 1717 | } else if (method == ClpSolve::tryDantzigWolfe) { |
| 1718 | abort(); |
| 1719 | } |
| 1720 | if (method == ClpSolve::usePrimalorSprint) { |
| 1721 | if (doSprint < 0) { |
| 1722 | if (numberElements < 500000) { |
| 1723 | // Small problem |
| 1724 | if (numberRows * 10 > numberColumns || numberColumns < 6000 |
| 1725 | || (numberRows * 20 > numberColumns && !plusMinus)) |
| 1726 | method = ClpSolve::usePrimal; // switch off sprint |
| 1727 | } else { |
| 1728 | // larger problem |
| 1729 | if (numberRows * 8 > numberColumns) |
| 1730 | method = ClpSolve::usePrimal; // switch off sprint |
| 1731 | // but make lightweight |
| 1732 | if (numberRows * 10 > numberColumns || numberColumns < 6000 |
| 1733 | || (numberRows * 20 > numberColumns && !plusMinus)) |
| 1734 | maxSprintPass = 10; |
| 1735 | } |
| 1736 | } else if (doSprint == 0) { |
| 1737 | method = ClpSolve::usePrimal; // switch off sprint |
| 1738 | } |
| 1739 | } |
| 1740 | if (method == ClpSolve::useDual) { |
| 1741 | #ifdef CLP_USEFUL_PRINTOUT |
| 1742 | debugInt[6] = 1; |
| 1743 | #endif |
| 1744 | double *saveLower = NULL; |
| 1745 | double *saveUpper = NULL; |
| 1746 | if (presolve == ClpSolve::presolveOn) { |
| 1747 | int numberInfeasibilities = model2->tightenPrimalBounds(0.0, 0); |
| 1748 | if (numberInfeasibilities) { |
| 1749 | handler_->message(CLP_INFEASIBLE, messages_) |
| 1750 | << CoinMessageEol; |
| 1751 | delete model2; |
| 1752 | model2 = this; |
| 1753 | presolve = ClpSolve::presolveOff; |
| 1754 | } |
| 1755 | } else if (numberRows_ + numberColumns_ > 5000) { |
| 1756 | // do anyway |
| 1757 | saveLower = new double[numberRows_ + numberColumns_]; |
| 1758 | CoinMemcpyN(model2->columnLower(), numberColumns_, saveLower); |
| 1759 | CoinMemcpyN(model2->rowLower(), numberRows_, saveLower + numberColumns_); |
| 1760 | saveUpper = new double[numberRows_ + numberColumns_]; |
| 1761 | CoinMemcpyN(model2->columnUpper(), numberColumns_, saveUpper); |
| 1762 | CoinMemcpyN(model2->rowUpper(), numberRows_, saveUpper + numberColumns_); |
| 1763 | int numberInfeasibilities = model2->tightenPrimalBounds(); |
| 1764 | if (numberInfeasibilities) { |
| 1765 | handler_->message(CLP_INFEASIBLE, messages_) |
| 1766 | << CoinMessageEol; |
| 1767 | CoinMemcpyN(saveLower, numberColumns_, model2->columnLower()); |
| 1768 | CoinMemcpyN(saveLower + numberColumns_, numberRows_, model2->rowLower()); |
| 1769 | delete[] saveLower; |
| 1770 | saveLower = NULL; |
| 1771 | CoinMemcpyN(saveUpper, numberColumns_, model2->columnUpper()); |
| 1772 | CoinMemcpyN(saveUpper + numberColumns_, numberRows_, model2->rowUpper()); |
| 1773 | delete[] saveUpper; |
| 1774 | saveUpper = NULL; |
| 1775 | // return if wanted |
| 1776 | if (options.infeasibleReturn() || (moreSpecialOptions_ & 1) != 0) |
| 1777 | return -1; |
| 1778 | } |
| 1779 | } |
| 1780 | #ifndef COIN_HAS_VOL |
| 1781 | // switch off idiot and volume for now |
| 1782 | doIdiot = 0; |
| 1783 | #endif |
| 1784 | // pick up number passes |
| 1785 | int nPasses = 0; |
| 1786 | int numberNotE = 0; |
| 1787 | #ifndef SLIM_CLP |
| 1788 | if ((doIdiot < 0 && plusMinus) || doIdiot > 0) { |
| 1789 | // See if candidate for idiot |
| 1790 | nPasses = 0; |
| 1791 | Idiot info(*model2); |
| 1792 | info.setStrategy(idiotOptions | info.getStrategy()); |
| 1793 | // Get average number of elements per column |
| 1794 | double ratio = static_cast< double >(numberElements) / static_cast< double >(numberColumns); |
| 1795 | // look at rhs |
| 1796 | int iRow; |
| 1797 | double largest = 0.0; |
| 1798 | double smallest = 1.0e30; |
| 1799 | for (iRow = 0; iRow < numberRows; iRow++) { |
| 1800 | double value1 = model2->rowLower_[iRow]; |
| 1801 | if (value1 && value1 > -1.0e31) { |
| 1802 | largest = CoinMax(largest, fabs(value1)); |
| 1803 | smallest = CoinMin(smallest, fabs(value1)); |
| 1804 | } |
| 1805 | double value2 = model2->rowUpper_[iRow]; |
| 1806 | if (value2 && value2 < 1.0e31) { |
| 1807 | largest = CoinMax(largest, fabs(value2)); |
| 1808 | smallest = CoinMin(smallest, fabs(value2)); |
| 1809 | } |
| 1810 | if (value2 > value1) { |
| 1811 | numberNotE++; |
| 1812 | } |
| 1813 | } |
| 1814 | if (doIdiot < 0) { |
| 1815 | if (numberRows > 200 && numberColumns > 5000 && ratio >= 3.0 && largest / smallest < 1.1 && !numberNotE) { |
| 1816 | nPasses = 71; |
| 1817 | } |
| 1818 | } |
| 1819 | if (doIdiot > 0) { |
| 1820 | nPasses = CoinMax(nPasses, doIdiot); |
| 1821 | if (nPasses > 70) { |
| 1822 | info.setStartingWeight(1.0e3); |
| 1823 | info.setDropEnoughFeasibility(0.01); |
| 1824 | } |
| 1825 | } |
| 1826 | if (nPasses > 20) { |
| 1827 | #ifdef COIN_HAS_VOL |
| 1828 | int returnCode = solveWithVolume(model2, nPasses, saveDoIdiot); |
| 1829 | if (!returnCode) { |
| 1830 | time2 = CoinCpuTime(); |
| 1831 | timeIdiot = time2 - timeX; |
| 1832 | handler_->message(CLP_INTERVAL_TIMING, messages_) |
| 1833 | << "Idiot Crash" << timeIdiot << time2 - time1 |
| 1834 | << CoinMessageEol; |
| 1835 | timeX = time2; |
| 1836 | } else { |
| 1837 | nPasses = 0; |
| 1838 | } |
1474 | | if (matrix_->type() == 11) { |
1475 | | // network - switch off stuff |
1476 | | doIdiot = 0; |
1477 | | //doSprint=0; |
1478 | | } |
1479 | | #endif |
1480 | | #endif |
1481 | | int numberColumns = model2->numberColumns(); |
1482 | | int numberRows = model2->numberRows(); |
1483 | | // If not all slack basis - switch off all except sprint |
1484 | | int numberRowsBasic = 0; |
1485 | | int iRow; |
1486 | | for (iRow = 0; iRow < numberRows; iRow++) |
1487 | | if (model2->getRowStatus(iRow) == basic) |
1488 | | numberRowsBasic++; |
1489 | | if (numberRowsBasic < numberRows && objective_->type()<2) { |
1490 | | doIdiot = 0; |
1491 | | doCrash = 0; |
1492 | | //doSprint=0; |
1493 | | } |
1494 | | if (options.getSpecialOption(3) == 0) { |
1495 | | if(numberElements > 100000) |
1496 | | plusMinus = true; |
1497 | | if(numberElements > 10000 && (doIdiot || doSprint)) |
1498 | | plusMinus = true; |
1499 | | } else if ((specialOptions_ & 1024) != 0) { |
1500 | | plusMinus = true; |
1501 | | } |
| 1840 | nPasses = 0; |
| 1841 | #endif |
| 1842 | } else { |
| 1843 | nPasses = 0; |
| 1844 | } |
| 1845 | } |
| 1846 | #endif |
| 1847 | if (doCrash) { |
| 1848 | #ifdef ABC_INHERIT |
| 1849 | if (!model2->abcState()) { |
| 1850 | #endif |
| 1851 | switch (doCrash) { |
| 1852 | // standard |
| 1853 | case 1: |
| 1854 | model2->crash(1000, 1); |
| 1855 | break; |
| 1856 | // As in paper by Solow and Halim (approx) |
| 1857 | case 2: |
| 1858 | case 3: |
| 1859 | model2->crash(model2->dualBound(), 0); |
| 1860 | break; |
| 1861 | // Just put free in basis |
| 1862 | case 4: |
| 1863 | model2->crash(0.0, 3); |
| 1864 | break; |
| 1865 | } |
| 1866 | #ifdef ABC_INHERIT |
| 1867 | } else if (doCrash >= 0) { |
| 1868 | model2->setAbcState(model2->abcState() | 256 * doCrash); |
| 1869 | } |
| 1870 | #endif |
| 1871 | } |
| 1872 | if (!nPasses) { |
| 1873 | int saveOptions = model2->specialOptions(); |
| 1874 | if (model2->numberRows() > 100) |
| 1875 | model2->setSpecialOptions(saveOptions | 64); // go as far as possible |
| 1876 | //int numberRows = model2->numberRows(); |
| 1877 | //int numberColumns = model2->numberColumns(); |
| 1878 | if (dynamic_cast< ClpPackedMatrix * >(matrix_)) { |
| 1879 | // See if original wanted vector |
| 1880 | ClpPackedMatrix *clpMatrixO = dynamic_cast< ClpPackedMatrix * >(matrix_); |
| 1881 | ClpMatrixBase *matrix = model2->clpMatrix(); |
| 1882 | if (dynamic_cast< ClpPackedMatrix * >(matrix) && clpMatrixO->wantsSpecialColumnCopy()) { |
| 1883 | ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(matrix); |
| 1884 | clpMatrix->makeSpecialColumnCopy(); |
| 1885 | //model2->setSpecialOptions(model2->specialOptions()|256); // to say no row copy for comparisons |
| 1886 | model2->dual(0); |
| 1887 | clpMatrix->releaseSpecialColumnCopy(); |
| 1888 | } else { |
| 1889 | #ifndef ABC_INHERIT |
| 1890 | model2->dual(0); |
| 1891 | #else |
| 1892 | model2->dealWithAbc(0, 0, interrupt); |
| 1893 | #endif |
| 1894 | } |
| 1895 | } else { |
| 1896 | model2->dual(0); |
| 1897 | } |
| 1898 | } else if (!numberNotE && 0) { |
| 1899 | // E so we can do in another way |
| 1900 | double *pi = model2->dualRowSolution(); |
| 1901 | int i; |
| 1902 | int numberColumns = model2->numberColumns(); |
| 1903 | int numberRows = model2->numberRows(); |
| 1904 | double *saveObj = new double[numberColumns]; |
| 1905 | CoinMemcpyN(model2->objective(), numberColumns, saveObj); |
| 1906 | CoinMemcpyN(model2->objective(), |
| 1907 | numberColumns, model2->dualColumnSolution()); |
| 1908 | model2->clpMatrix()->transposeTimes(-1.0, pi, model2->dualColumnSolution()); |
| 1909 | CoinMemcpyN(model2->dualColumnSolution(), |
| 1910 | numberColumns, model2->objective()); |
| 1911 | const double *rowsol = model2->primalRowSolution(); |
| 1912 | double offset = 0.0; |
| 1913 | for (i = 0; i < numberRows; i++) { |
| 1914 | offset += pi[i] * rowsol[i]; |
| 1915 | } |
| 1916 | double value2; |
| 1917 | model2->getDblParam(ClpObjOffset, value2); |
| 1918 | //printf("Offset %g %g\n",offset,value2); |
| 1919 | model2->setDblParam(ClpObjOffset, value2 - offset); |
| 1920 | model2->setPerturbation(51); |
| 1921 | //model2->setRowObjective(pi); |
| 1922 | // zero out pi |
| 1923 | //memset(pi,0,numberRows*sizeof(double)); |
| 1924 | // Could put some in basis - only partially tested |
| 1925 | model2->allSlackBasis(); |
| 1926 | //model2->factorization()->maximumPivots(200); |
| 1927 | //model2->setLogLevel(63); |
| 1928 | // solve |
| 1929 | model2->dual(0); |
| 1930 | model2->setDblParam(ClpObjOffset, value2); |
| 1931 | CoinMemcpyN(saveObj, numberColumns, model2->objective()); |
| 1932 | // zero out pi |
| 1933 | //memset(pi,0,numberRows*sizeof(double)); |
| 1934 | //model2->setRowObjective(pi); |
| 1935 | delete[] saveObj; |
| 1936 | //model2->dual(0); |
| 1937 | model2->setPerturbation(50); |
| 1938 | model2->primal(); |
| 1939 | } else { |
| 1940 | // solve |
| 1941 | model2->setPerturbation(100); |
| 1942 | model2->dual(2); |
| 1943 | model2->setPerturbation(50); |
| 1944 | model2->dual(0); |
| 1945 | } |
| 1946 | if (saveLower) { |
| 1947 | CoinMemcpyN(saveLower, numberColumns_, model2->columnLower()); |
| 1948 | CoinMemcpyN(saveLower + numberColumns_, numberRows_, model2->rowLower()); |
| 1949 | delete[] saveLower; |
| 1950 | saveLower = NULL; |
| 1951 | CoinMemcpyN(saveUpper, numberColumns_, model2->columnUpper()); |
| 1952 | CoinMemcpyN(saveUpper + numberColumns_, numberRows_, model2->rowUpper()); |
| 1953 | delete[] saveUpper; |
| 1954 | saveUpper = NULL; |
| 1955 | } |
| 1956 | time2 = CoinCpuTime(); |
| 1957 | timeCore = time2 - timeX; |
| 1958 | handler_->message(CLP_INTERVAL_TIMING, messages_) |
| 1959 | << "Dual" << timeCore << time2 - time1 |
| 1960 | << CoinMessageEol; |
| 1961 | timeX = time2; |
| 1962 | } else if (method == ClpSolve::usePrimal) { |
| 1963 | #ifdef CLP_USEFUL_PRINTOUT |
| 1964 | debugInt[6] = 2; |
| 1965 | #endif |
1503 | | // Statistics (+1,-1, other) - used to decide on strategy if not +-1 |
1504 | | CoinBigIndex statistics[3] = { -1, 0, 0}; |
1505 | | if (plusMinus) { |
1506 | | saveMatrix = model2->clpMatrix(); |
1507 | | #ifndef NO_RTTI |
1508 | | ClpPackedMatrix* clpMatrix = |
1509 | | dynamic_cast< ClpPackedMatrix*>(saveMatrix); |
1510 | | #else |
1511 | | ClpPackedMatrix* clpMatrix = NULL; |
1512 | | if (saveMatrix->type() == 1) |
1513 | | clpMatrix = |
1514 | | static_cast< ClpPackedMatrix*>(saveMatrix); |
1515 | | #endif |
1516 | | if (clpMatrix) { |
1517 | | ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix())); |
1518 | | if (newMatrix->getIndices()) { |
1519 | | // CHECKME This test of specialOptions and the one above |
1520 | | // don't seem compatible. |
1521 | | #ifndef ABC_INHERIT |
1522 | | if ((specialOptions_ & 1024) == 0) { |
1523 | | model2->replaceMatrix(newMatrix); |
1524 | | } else { |
1525 | | #endif |
1526 | | // in integer (or abc) - just use for sprint/idiot |
1527 | | saveMatrix = NULL; |
1528 | | delete newMatrix; |
1529 | | #ifndef ABC_INHERIT |
1530 | | } |
1531 | | #endif |
1532 | | } else { |
1533 | | handler_->message(CLP_MATRIX_CHANGE, messages_) |
1534 | | << "+- 1" |
1535 | | << CoinMessageEol; |
1536 | | CoinMemcpyN(newMatrix->startPositive(), 3, statistics); |
1537 | | saveMatrix = NULL; |
1538 | | plusMinus = false; |
1539 | | delete newMatrix; |
1540 | | } |
| 1967 | if (doIdiot) { |
| 1968 | int nPasses = 0; |
| 1969 | Idiot info(*model2); |
| 1970 | info.setStrategy(idiotOptions | info.getStrategy()); |
| 1971 | // Get average number of elements per column |
| 1972 | double ratio = static_cast< double >(numberElements) / static_cast< double >(numberColumns); |
| 1973 | // look at rhs |
| 1974 | int iRow; |
| 1975 | double largest = 0.0; |
| 1976 | double smallest = 1.0e30; |
| 1977 | double largestGap = 0.0; |
| 1978 | int numberNotE = 0; |
| 1979 | for (iRow = 0; iRow < numberRows; iRow++) { |
| 1980 | double value1 = model2->rowLower_[iRow]; |
| 1981 | if (value1 && value1 > -1.0e31) { |
| 1982 | largest = CoinMax(largest, fabs(value1)); |
| 1983 | smallest = CoinMin(smallest, fabs(value1)); |
| 1984 | } |
| 1985 | double value2 = model2->rowUpper_[iRow]; |
| 1986 | if (value2 && value2 < 1.0e31) { |
| 1987 | largest = CoinMax(largest, fabs(value2)); |
| 1988 | smallest = CoinMin(smallest, fabs(value2)); |
| 1989 | } |
| 1990 | if (value2 > value1) { |
| 1991 | numberNotE++; |
| 1992 | if (value2 > 1.0e31 || value1 < -1.0e31) |
| 1993 | largestGap = COIN_DBL_MAX; |
| 1994 | else |
| 1995 | largestGap = value2 - value1; |
| 1996 | } |
| 1997 | } |
| 1998 | bool increaseSprint = plusMinus; |
| 1999 | if ((specialOptions_ & 1024) != 0) |
| 2000 | increaseSprint = false; |
| 2001 | if (!plusMinus) { |
| 2002 | // If 90% +- 1 then go for sprint |
| 2003 | if (statistics[0] >= 0 && 10 * statistics[2] < statistics[0] + statistics[1]) |
| 2004 | increaseSprint = true; |
| 2005 | } |
| 2006 | int tryIt = tryItSave; |
| 2007 | if (numberRows < 1000 && numberColumns < 3000) |
| 2008 | tryIt = 0; |
| 2009 | if (tryIt) { |
| 2010 | if (increaseSprint) { |
| 2011 | info.setStartingWeight(1.0e3); |
| 2012 | info.setReduceIterations(6); |
| 2013 | // also be more lenient on infeasibilities |
| 2014 | info.setDropEnoughFeasibility(0.5 * info.getDropEnoughFeasibility()); |
| 2015 | info.setDropEnoughWeighted(-2.0); |
| 2016 | if (largest / smallest > 2.0) { |
| 2017 | nPasses = 10 + numberColumns / 100000; |
| 2018 | nPasses = CoinMin(nPasses, 50); |
| 2019 | nPasses = CoinMax(nPasses, 15); |
| 2020 | if (numberRows > 20000 && nPasses > 5) { |
| 2021 | // Might as well go for it |
| 2022 | nPasses = CoinMax(nPasses, 71); |
| 2023 | } else if (numberRows > 2000 && nPasses > 5) { |
| 2024 | nPasses = CoinMax(nPasses, 50); |
| 2025 | } else if (numberElements < 3 * numberColumns) { |
| 2026 | nPasses = CoinMin(nPasses, 10); // probably not worh it |
| 2027 | if (doIdiot < 0) |
| 2028 | info.setLightweight(1); // say lightweight idiot |
| 2029 | } else { |
| 2030 | if (doIdiot < 0) |
| 2031 | info.setLightweight(1); // say lightweight idiot |
| 2032 | } |
| 2033 | } else if (largest / smallest > 1.01 || numberElements <= 3 * numberColumns) { |
| 2034 | nPasses = 10 + numberColumns / 1000; |
| 2035 | nPasses = CoinMin(nPasses, 100); |
| 2036 | nPasses = CoinMax(nPasses, 30); |
| 2037 | if (numberRows > 25000) { |
| 2038 | // Might as well go for it |
| 2039 | nPasses = CoinMax(nPasses, 71); |
| 2040 | } |
| 2041 | if (!largestGap) |
| 2042 | nPasses *= 2; |
1559 | | // only do primal if sprint or idiot |
1560 | | if (doSprint > 0) { |
1561 | | method = ClpSolve::usePrimalorSprint; |
1562 | | } else if (doIdiot > 0) { |
1563 | | method = ClpSolve::usePrimal; |
1564 | | } else { |
1565 | | if (numberElements < 500000) { |
1566 | | // Small problem |
1567 | | if(numberRows * 10 > numberColumns || numberColumns < 6000 |
1568 | | || (numberRows * 20 > numberColumns && !plusMinus)) |
1569 | | doSprint = 0; // switch off sprint |
1570 | | } else { |
1571 | | // larger problem |
1572 | | if(numberRows * 8 > numberColumns) |
1573 | | doSprint = 0; // switch off sprint |
1574 | | } |
1575 | | // switch off idiot or sprint if any free variable |
1576 | | // switch off sprint if very few with costs |
1577 | | // or great variation in cost |
1578 | | int iColumn; |
1579 | | const double * columnLower = model2->columnLower(); |
1580 | | const double * columnUpper = model2->columnUpper(); |
1581 | | const double * objective = model2->objective(); |
1582 | | int nObj=0; |
1583 | | int nFree=0; |
1584 | | double smallestObj=COIN_DBL_MAX; |
1585 | | double largestObj=0.0; |
1586 | | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
1587 | | if (columnLower[iColumn] < -1.0e10 && columnUpper[iColumn] > 1.0e10) { |
1588 | | nFree++; |
1589 | | } else if (objective[iColumn]) { |
1590 | | nObj++; |
1591 | | smallestObj=CoinMin(smallestObj,objective[iColumn]); |
1592 | | largestObj=CoinMax(largestObj,objective[iColumn]); |
1593 | | } |
1594 | | } |
1595 | | if (nObj*10<numberColumns || |
1596 | | smallestObj*10.0<largestObj) |
1597 | | doSprint=0; |
1598 | | if (nFree) |
1599 | | doIdiot=0; |
1600 | | if (nFree*10>numberRows) |
1601 | | doSprint=0; |
1602 | | int nPasses = 0; |
1603 | | // look at rhs |
1604 | | int iRow; |
1605 | | double largest = 0.0; |
1606 | | double smallest = 1.0e30; |
1607 | | double largestGap = 0.0; |
1608 | | int numberNotE = 0; |
1609 | | bool notInteger = false; |
1610 | | for (iRow = 0; iRow < numberRows; iRow++) { |
1611 | | double value1 = model2->rowLower_[iRow]; |
1612 | | if (value1 && value1 > -1.0e31) { |
1613 | | largest = CoinMax(largest, fabs(value1)); |
1614 | | smallest = CoinMin(smallest, fabs(value1)); |
1615 | | if (fabs(value1 - floor(value1 + 0.5)) > 1.0e-8) { |
1616 | | notInteger = true; |
1617 | | break; |
1618 | | } |
1619 | | } |
1620 | | double value2 = model2->rowUpper_[iRow]; |
1621 | | if (value2 && value2 < 1.0e31) { |
1622 | | largest = CoinMax(largest, fabs(value2)); |
1623 | | smallest = CoinMin(smallest, fabs(value2)); |
1624 | | if (fabs(value2 - floor(value2 + 0.5)) > 1.0e-8) { |
1625 | | notInteger = true; |
1626 | | break; |
1627 | | } |
1628 | | } |
1629 | | // CHECKME This next bit can't be right... |
1630 | | if (value2 > value1) { |
1631 | | numberNotE++; |
1632 | | //if (value2 > 1.0e31 || value1 < -1.0e31) |
1633 | | // largestGap = COIN_DBL_MAX; |
1634 | | //else |
1635 | | // largestGap = value2 - value1; |
1636 | | } |
1637 | | } |
1638 | | int tryIt = (numberRows > 200 && numberColumns > 2000 && |
1639 | | (numberColumns > 2 * numberRows || (method != ClpSolve::useDual && (specialOptions_ & 1024) != 0))) ? 3 : 0; |
1640 | | tryItSave = tryIt; |
1641 | | if (numberRows < 1000 && numberColumns < 3000) |
1642 | | tryIt = 0; |
1643 | | if (notInteger) |
1644 | | tryIt = 0; |
1645 | | if (largest / smallest > 10 || (largest / smallest > 2.0 && largest > 50)) |
1646 | | tryIt = 0; |
1647 | | if (tryIt) { |
1648 | | if (largest / smallest > 2.0) { |
1649 | | nPasses = 10 + numberColumns / 100000; |
1650 | | nPasses = CoinMin(nPasses, 50); |
1651 | | nPasses = CoinMax(nPasses, 15); |
1652 | | if (numberRows > 20000 && nPasses > 5) { |
1653 | | // Might as well go for it |
1654 | | nPasses = CoinMax(nPasses, 71); |
1655 | | } else if (numberRows > 2000 && nPasses > 5) { |
1656 | | nPasses = CoinMax(nPasses, 50); |
1657 | | } else if (numberElements < 3 * numberColumns) { |
1658 | | nPasses = CoinMin(nPasses, 10); // probably not worh it |
1659 | | } |
1660 | | } else if (largest / smallest > 1.01 || numberElements <= 3 * numberColumns) { |
1661 | | nPasses = 10 + numberColumns / 1000; |
1662 | | nPasses = CoinMin(nPasses, 100); |
1663 | | nPasses = CoinMax(nPasses, 30); |
1664 | | if (numberRows > 25000) { |
1665 | | // Might as well go for it |
1666 | | nPasses = CoinMax(nPasses, 71); |
1667 | | } |
1668 | | if (!largestGap) |
1669 | | nPasses *= 2; |
1670 | | } else { |
1671 | | nPasses = 10 + numberColumns / 1000; |
1672 | | nPasses = CoinMax(nPasses, 100); |
1673 | | if (!largestGap) |
1674 | | nPasses *= 2; |
1675 | | nPasses = CoinMin(nPasses, 200); |
1676 | | } |
1677 | | } |
1678 | | //printf("%d rows %d cols plus %c tryIt %c largest %g smallest %g largestGap %g npasses %d sprint %c\n", |
1679 | | // numberRows,numberColumns,plusMinus ? 'Y' : 'N', |
1680 | | // tryIt ? 'Y' :'N',largest,smallest,largestGap,nPasses,doSprint ? 'Y' :'N'); |
1681 | | //exit(0); |
1682 | | if (!tryIt || nPasses <= 5) |
1683 | | doIdiot = 0; |
1684 | | #if CLPSOLVE_ACTIONS |
1685 | | if (doIdiot&&doSprint<0&&wasAutomatic && |
1686 | | 20*model2->numberRows()>model2->numberColumns()) |
1687 | | doSprint=0; // switch off sprint |
1688 | | #endif |
1689 | | if (doSprint) { |
1690 | | method = ClpSolve::usePrimalorSprint; |
1691 | | } else if (doIdiot) { |
1692 | | method = ClpSolve::usePrimal; |
1693 | | } else { |
1694 | | method = ClpSolve::useDual; |
1695 | | } |
1696 | | } |
| 2068 | if (numberRows > 25000 || numberColumns > 5 * numberRows) { |
| 2069 | nPasses = 50; |
| 2070 | info.setLightweight(0); // say not lightweight idiot |
| 2071 | } else if (numberColumns > 4 * numberRows) { |
| 2072 | nPasses = 20; |
| 2073 | } else { |
| 2074 | nPasses = 15; |
| 2075 | } |
1752 | | debugInt[6]=1; |
1753 | | #endif |
1754 | | double * saveLower = NULL; |
1755 | | double * saveUpper = NULL; |
1756 | | if (presolve == ClpSolve::presolveOn) { |
1757 | | int numberInfeasibilities = model2->tightenPrimalBounds(0.0, 0); |
1758 | | if (numberInfeasibilities) { |
1759 | | handler_->message(CLP_INFEASIBLE, messages_) |
1760 | | << CoinMessageEol; |
1761 | | delete model2; |
1762 | | model2 = this; |
1763 | | presolve = ClpSolve::presolveOff; |
1764 | | } |
1765 | | } else if (numberRows_ + numberColumns_ > 5000) { |
1766 | | // do anyway |
1767 | | saveLower = new double[numberRows_+numberColumns_]; |
1768 | | CoinMemcpyN(model2->columnLower(), numberColumns_, saveLower); |
1769 | | CoinMemcpyN(model2->rowLower(), numberRows_, saveLower + numberColumns_); |
1770 | | saveUpper = new double[numberRows_+numberColumns_]; |
1771 | | CoinMemcpyN(model2->columnUpper(), numberColumns_, saveUpper); |
1772 | | CoinMemcpyN(model2->rowUpper(), numberRows_, saveUpper + numberColumns_); |
1773 | | int numberInfeasibilities = model2->tightenPrimalBounds(); |
1774 | | if (numberInfeasibilities) { |
1775 | | handler_->message(CLP_INFEASIBLE, messages_) |
1776 | | << CoinMessageEol; |
1777 | | CoinMemcpyN(saveLower, numberColumns_, model2->columnLower()); |
1778 | | CoinMemcpyN(saveLower + numberColumns_, numberRows_, model2->rowLower()); |
1779 | | delete [] saveLower; |
1780 | | saveLower = NULL; |
1781 | | CoinMemcpyN(saveUpper, numberColumns_, model2->columnUpper()); |
1782 | | CoinMemcpyN(saveUpper + numberColumns_, numberRows_, model2->rowUpper()); |
1783 | | delete [] saveUpper; |
1784 | | saveUpper = NULL; |
1785 | | // return if wanted |
1786 | | if (options.infeasibleReturn() || |
1787 | | (moreSpecialOptions_ & 1) != 0) |
1788 | | return -1; |
1789 | | } |
| 2118 | debugInt[6] = 4; |
| 2119 | debugInt[7] = nPasses; |
| 2120 | #endif |
| 2121 | if (nPasses > 70) { |
| 2122 | info.setStartingWeight(1.0e3); |
| 2123 | info.setReduceIterations(6); |
| 2124 | //if (nPasses > 200) |
| 2125 | //info.setFeasibilityTolerance(1.0e-9); |
| 2126 | //if (nPasses > 1900) |
| 2127 | //info.setWeightFactor(0.93); |
| 2128 | if (nPasses > 900) { |
| 2129 | double reductions = nPasses / 6.0; |
| 2130 | if (nPasses < 5000) { |
| 2131 | reductions /= 12.0; |
| 2132 | } else { |
| 2133 | reductions /= 13.0; |
| 2134 | info.setStartingWeight(1.0e4); |
| 2135 | } |
| 2136 | double ratio = 1.0 / std::pow(10.0, (1.0 / reductions)); |
| 2137 | printf("%d passes reduction factor %g\n", nPasses, ratio); |
| 2138 | info.setWeightFactor(ratio); |
| 2139 | } else if (nPasses > 500) { |
| 2140 | info.setWeightFactor(0.7); |
| 2141 | } else if (nPasses > 200) { |
| 2142 | info.setWeightFactor(0.5); |
1884 | | if (!nPasses) { |
1885 | | int saveOptions = model2->specialOptions(); |
1886 | | if (model2->numberRows() > 100) |
1887 | | model2->setSpecialOptions(saveOptions | 64); // go as far as possible |
1888 | | //int numberRows = model2->numberRows(); |
1889 | | //int numberColumns = model2->numberColumns(); |
1890 | | if (dynamic_cast< ClpPackedMatrix*>(matrix_)) { |
1891 | | // See if original wanted vector |
1892 | | ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_); |
1893 | | ClpMatrixBase * matrix = model2->clpMatrix(); |
1894 | | if (dynamic_cast< ClpPackedMatrix*>(matrix) && clpMatrixO->wantsSpecialColumnCopy()) { |
1895 | | ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix); |
1896 | | clpMatrix->makeSpecialColumnCopy(); |
1897 | | //model2->setSpecialOptions(model2->specialOptions()|256); // to say no row copy for comparisons |
1898 | | model2->dual(0); |
1899 | | clpMatrix->releaseSpecialColumnCopy(); |
1900 | | } else { |
1901 | | #ifndef ABC_INHERIT |
1902 | | model2->dual(0); |
1903 | | #else |
1904 | | model2->dealWithAbc(0,0,interrupt); |
1905 | | #endif |
1906 | | } |
1907 | | } else { |
1908 | | model2->dual(0); |
1909 | | } |
1910 | | } else if (!numberNotE && 0) { |
1911 | | // E so we can do in another way |
1912 | | double * pi = model2->dualRowSolution(); |
1913 | | int i; |
1914 | | int numberColumns = model2->numberColumns(); |
1915 | | int numberRows = model2->numberRows(); |
1916 | | double * saveObj = new double[numberColumns]; |
1917 | | CoinMemcpyN(model2->objective(), numberColumns, saveObj); |
1918 | | CoinMemcpyN(model2->objective(), |
1919 | | numberColumns, model2->dualColumnSolution()); |
1920 | | model2->clpMatrix()->transposeTimes(-1.0, pi, model2->dualColumnSolution()); |
1921 | | CoinMemcpyN(model2->dualColumnSolution(), |
1922 | | numberColumns, model2->objective()); |
1923 | | const double * rowsol = model2->primalRowSolution(); |
1924 | | double offset = 0.0; |
1925 | | for (i = 0; i < numberRows; i++) { |
1926 | | offset += pi[i] * rowsol[i]; |
1927 | | } |
1928 | | double value2; |
1929 | | model2->getDblParam(ClpObjOffset, value2); |
1930 | | //printf("Offset %g %g\n",offset,value2); |
1931 | | model2->setDblParam(ClpObjOffset, value2 - offset); |
1932 | | model2->setPerturbation(51); |
1933 | | //model2->setRowObjective(pi); |
1934 | | // zero out pi |
1935 | | //memset(pi,0,numberRows*sizeof(double)); |
1936 | | // Could put some in basis - only partially tested |
1937 | | model2->allSlackBasis(); |
1938 | | //model2->factorization()->maximumPivots(200); |
1939 | | //model2->setLogLevel(63); |
1940 | | // solve |
1941 | | model2->dual(0); |
1942 | | model2->setDblParam(ClpObjOffset, value2); |
1943 | | CoinMemcpyN(saveObj, numberColumns, model2->objective()); |
1944 | | // zero out pi |
1945 | | //memset(pi,0,numberRows*sizeof(double)); |
1946 | | //model2->setRowObjective(pi); |
1947 | | delete [] saveObj; |
1948 | | //model2->dual(0); |
1949 | | model2->setPerturbation(50); |
1950 | | model2->primal(); |
1951 | | } else { |
1952 | | // solve |
1953 | | model2->setPerturbation(100); |
1954 | | model2->dual(2); |
1955 | | model2->setPerturbation(50); |
1956 | | model2->dual(0); |
1957 | | } |
1958 | | if (saveLower) { |
1959 | | CoinMemcpyN(saveLower, numberColumns_, model2->columnLower()); |
1960 | | CoinMemcpyN(saveLower + numberColumns_, numberRows_, model2->rowLower()); |
1961 | | delete [] saveLower; |
1962 | | saveLower = NULL; |
1963 | | CoinMemcpyN(saveUpper, numberColumns_, model2->columnUpper()); |
1964 | | CoinMemcpyN(saveUpper + numberColumns_, numberRows_, model2->rowUpper()); |
1965 | | delete [] saveUpper; |
1966 | | saveUpper = NULL; |
1967 | | } |
1968 | | time2 = CoinCpuTime(); |
1969 | | timeCore = time2 - timeX; |
1970 | | handler_->message(CLP_INTERVAL_TIMING, messages_) |
1971 | | << "Dual" << timeCore << time2 - time1 |
1972 | | << CoinMessageEol; |
1973 | | timeX = time2; |
1974 | | } else if (method == ClpSolve::usePrimal) { |
1975 | | #ifdef CLP_USEFUL_PRINTOUT |
1976 | | debugInt[6]=2; |
1977 | | #endif |
1978 | | #ifndef SLIM_CLP |
1979 | | if (doIdiot) { |
1980 | | int nPasses = 0; |
1981 | | Idiot info(*model2); |
1982 | | info.setStrategy(idiotOptions | info.getStrategy()); |
1983 | | // Get average number of elements per column |
1984 | | double ratio = static_cast<double> (numberElements) / static_cast<double> (numberColumns); |
1985 | | // look at rhs |
1986 | | int iRow; |
1987 | | double largest = 0.0; |
1988 | | double smallest = 1.0e30; |
1989 | | double largestGap = 0.0; |
1990 | | int numberNotE = 0; |
1991 | | for (iRow = 0; iRow < numberRows; iRow++) { |
1992 | | double value1 = model2->rowLower_[iRow]; |
1993 | | if (value1 && value1 > -1.0e31) { |
1994 | | largest = CoinMax(largest, fabs(value1)); |
1995 | | smallest = CoinMin(smallest, fabs(value1)); |
1996 | | } |
1997 | | double value2 = model2->rowUpper_[iRow]; |
1998 | | if (value2 && value2 < 1.0e31) { |
1999 | | largest = CoinMax(largest, fabs(value2)); |
2000 | | smallest = CoinMin(smallest, fabs(value2)); |
2001 | | } |
2002 | | if (value2 > value1) { |
2003 | | numberNotE++; |
2004 | | if (value2 > 1.0e31 || value1 < -1.0e31) |
2005 | | largestGap = COIN_DBL_MAX; |
2006 | | else |
2007 | | largestGap = value2 - value1; |
2008 | | } |
2009 | | } |
2010 | | bool increaseSprint = plusMinus; |
2011 | | if ((specialOptions_ & 1024) != 0) |
2012 | | increaseSprint = false; |
2013 | | if (!plusMinus) { |
2014 | | // If 90% +- 1 then go for sprint |
2015 | | if (statistics[0] >= 0 && 10 * statistics[2] < statistics[0] + statistics[1]) |
2016 | | increaseSprint = true; |
2017 | | } |
2018 | | int tryIt = tryItSave; |
2019 | | if (numberRows < 1000 && numberColumns < 3000) |
2020 | | tryIt = 0; |
2021 | | if (tryIt) { |
2022 | | if (increaseSprint) { |
2023 | | info.setStartingWeight(1.0e3); |
2024 | | info.setReduceIterations(6); |
2025 | | // also be more lenient on infeasibilities |
2026 | | info.setDropEnoughFeasibility(0.5 * info.getDropEnoughFeasibility()); |
2027 | | info.setDropEnoughWeighted(-2.0); |
2028 | | if (largest / smallest > 2.0) { |
2029 | | nPasses = 10 + numberColumns / 100000; |
2030 | | nPasses = CoinMin(nPasses, 50); |
2031 | | nPasses = CoinMax(nPasses, 15); |
2032 | | if (numberRows > 20000 && nPasses > 5) { |
2033 | | // Might as well go for it |
2034 | | nPasses = CoinMax(nPasses, 71); |
2035 | | } else if (numberRows > 2000 && nPasses > 5) { |
2036 | | nPasses = CoinMax(nPasses, 50); |
2037 | | } else if (numberElements < 3 * numberColumns) { |
2038 | | nPasses = CoinMin(nPasses, 10); // probably not worh it |
2039 | | if (doIdiot < 0) |
2040 | | info.setLightweight(1); // say lightweight idiot |
2041 | | } else { |
2042 | | if (doIdiot < 0) |
2043 | | info.setLightweight(1); // say lightweight idiot |
2044 | | } |
2045 | | } else if (largest / smallest > 1.01 || numberElements <= 3 * numberColumns) { |
2046 | | nPasses = 10 + numberColumns / 1000; |
2047 | | nPasses = CoinMin(nPasses, 100); |
2048 | | nPasses = CoinMax(nPasses, 30); |
2049 | | if (numberRows > 25000) { |
2050 | | // Might as well go for it |
2051 | | nPasses = CoinMax(nPasses, 71); |
2052 | | } |
2053 | | if (!largestGap) |
2054 | | nPasses *= 2; |
2055 | | } else { |
2056 | | nPasses = 10 + numberColumns / 1000; |
2057 | | nPasses = CoinMin(nPasses, 200); |
2058 | | nPasses = CoinMax(nPasses, 100); |
2059 | | info.setStartingWeight(1.0e-1); |
2060 | | info.setReduceIterations(6); |
2061 | | if (!largestGap && nPasses <= 50) |
2062 | | nPasses *= 2; |
2063 | | //info.setFeasibilityTolerance(1.0e-7); |
2064 | | } |
2065 | | // If few passes - don't bother |
2066 | | if (nPasses <= 5 && !plusMinus) |
2067 | | nPasses = 0; |
2068 | | } else { |
2069 | | if (doIdiot < 0) |
2070 | | info.setLightweight(1); // say lightweight idiot |
2071 | | if (largest / smallest > 1.01 || numberNotE || statistics[2] > statistics[0] + statistics[1]) { |
2072 | | if (numberRows > 25000 || numberColumns > 5 * numberRows) { |
2073 | | nPasses = 50; |
2074 | | } else if (numberColumns > 4 * numberRows) { |
2075 | | nPasses = 20; |
2076 | | } else { |
2077 | | nPasses = 5; |
2078 | | } |
2079 | | } else { |
2080 | | if (numberRows > 25000 || numberColumns > 5 * numberRows) { |
2081 | | nPasses = 50; |
2082 | | info.setLightweight(0); // say not lightweight idiot |
2083 | | } else if (numberColumns > 4 * numberRows) { |
2084 | | nPasses = 20; |
2085 | | } else { |
2086 | | nPasses = 15; |
2087 | | } |
2088 | | } |
2089 | | if (ratio < 3.0) { |
2090 | | nPasses = static_cast<int> (ratio * static_cast<double> (nPasses) / 4.0); // probably not worth it |
2091 | | } else { |
2092 | | nPasses = CoinMax(nPasses, 5); |
2093 | | } |
2094 | | if (numberRows > 25000 && nPasses > 5) { |
2095 | | // Might as well go for it |
2096 | | nPasses = CoinMax(nPasses, 71); |
2097 | | } else if (increaseSprint) { |
2098 | | nPasses *= 2; |
2099 | | nPasses = CoinMin(nPasses, 71); |
2100 | | } else if (nPasses == 5 && ratio > 5.0) { |
2101 | | nPasses = static_cast<int> (static_cast<double> (nPasses) * (ratio / 5.0)); // increase if lots of elements per column |
2102 | | } |
2103 | | if (nPasses <= 5 && !plusMinus) |
2104 | | nPasses = 0; |
2105 | | //info.setStartingWeight(1.0e-1); |
2106 | | } |
2107 | | if (tryIt==1) { |
2108 | | idiotOptions |= 262144; |
2109 | | info.setStrategy(idiotOptions | info.getStrategy()); |
2110 | | //model2->setSpecialOptions(model2->specialOptions() |
2111 | | // |8388608); |
2112 | | } |
2113 | | } |
2114 | | if (doIdiot > 0) { |
2115 | | // pick up number passes |
2116 | | nPasses = options.getExtraInfo(1) % 1000000; |
2117 | | #ifdef COIN_HAS_VOL |
2118 | | int returnCode = solveWithVolume(model2, nPasses, saveDoIdiot); |
2119 | | nPasses=0; |
2120 | | if (!returnCode) { |
2121 | | time2 = CoinCpuTime(); |
2122 | | timeIdiot = time2 - timeX; |
2123 | | handler_->message(CLP_INTERVAL_TIMING, messages_) |
2124 | | << "Idiot Crash" << timeIdiot << time2 - time1 |
2125 | | << CoinMessageEol; |
2126 | | timeX = time2; |
2127 | | } |
2128 | | #endif |
2129 | | #ifdef CLP_USEFUL_PRINTOUT |
2130 | | debugInt[6]=4; |
2131 | | debugInt[7]=nPasses; |
2132 | | #endif |
2133 | | if (nPasses > 70) { |
2134 | | info.setStartingWeight(1.0e3); |
2135 | | info.setReduceIterations(6); |
2136 | | //if (nPasses > 200) |
2137 | | //info.setFeasibilityTolerance(1.0e-9); |
2138 | | //if (nPasses > 1900) |
2139 | | //info.setWeightFactor(0.93); |
2140 | | if (nPasses > 900) { |
2141 | | double reductions=nPasses/6.0; |
2142 | | if (nPasses<5000) { |
2143 | | reductions /= 12.0; |
2144 | | } else { |
2145 | | reductions /= 13.0; |
2146 | | info.setStartingWeight(1.0e4); |
2147 | | } |
2148 | | double ratio=1.0/std::pow(10.0,(1.0/reductions)); |
2149 | | printf("%d passes reduction factor %g\n",nPasses,ratio); |
2150 | | info.setWeightFactor(ratio); |
2151 | | } else if (nPasses > 500) { |
2152 | | info.setWeightFactor(0.7); |
2153 | | } else if (nPasses > 200) { |
2154 | | info.setWeightFactor(0.5); |
2155 | | } |
2156 | | if (maximumIterations()<nPasses) { |
2157 | | printf("Presuming maximumIterations is just for Idiot\n"); |
2158 | | nPasses=maximumIterations(); |
2159 | | setMaximumIterations(COIN_INT_MAX); |
2160 | | model2->setMaximumIterations(COIN_INT_MAX); |
2161 | | } |
2162 | | if (nPasses >= 10000&&nPasses<100000) { |
2163 | | int k = nPasses % 100; |
2164 | | nPasses /= 200; |
2165 | | info.setReduceIterations(3); |
2166 | | if (k) |
2167 | | info.setStartingWeight(1.0e2); |
2168 | | } |
2169 | | // also be more lenient on infeasibilities |
2170 | | info.setDropEnoughFeasibility(0.5 * info.getDropEnoughFeasibility()); |
2171 | | info.setDropEnoughWeighted(-2.0); |
2172 | | } else if (nPasses >= 50) { |
2173 | | info.setStartingWeight(1.0e3); |
2174 | | //info.setReduceIterations(6); |
2175 | | } |
2176 | | // For experimenting |
2177 | | if (nPasses < 70 && (nPasses % 10) > 0 && (nPasses % 10) < 4) { |
2178 | | info.setStartingWeight(1.0e3); |
2179 | | info.setLightweight(nPasses % 10); // special testing |
| 2157 | // also be more lenient on infeasibilities |
| 2158 | info.setDropEnoughFeasibility(0.5 * info.getDropEnoughFeasibility()); |
| 2159 | info.setDropEnoughWeighted(-2.0); |
| 2160 | } else if (nPasses >= 50) { |
| 2161 | info.setStartingWeight(1.0e3); |
| 2162 | //info.setReduceIterations(6); |
| 2163 | } |
| 2164 | // For experimenting |
| 2165 | if (nPasses < 70 && (nPasses % 10) > 0 && (nPasses % 10) < 4) { |
| 2166 | info.setStartingWeight(1.0e3); |
| 2167 | info.setLightweight(nPasses % 10); // special testing |
2259 | | // move duals and just copy primal |
2260 | | memcpy(temp.dualRowSolution(),dualModel2->primalColumnSolution(), |
2261 | | numberRows*sizeof(double)); |
2262 | | memcpy(temp.primalColumnSolution(),model2->primalColumnSolution(), |
2263 | | numberColumns*sizeof(double)); |
2264 | | #endif |
2265 | | delete dualModel2; |
2266 | | int numberRows=model2->numberRows(); |
2267 | | int numberColumns=model2->numberColumns(); |
2268 | | ClpSimplex * tempModel[2]; |
2269 | | tempModel[0]=model2; |
2270 | | tempModel[1]=&temp; |
2271 | | const double * primalColumn[2]; |
2272 | | const double * dualRow[2]; |
2273 | | double * dualColumn[2]; |
2274 | | double * primalRow[2]; |
2275 | | for (int i=0;i<2;i++) { |
2276 | | primalColumn[i]=tempModel[i]->primalColumnSolution(); |
2277 | | dualRow[i]=tempModel[i]->dualRowSolution(); |
2278 | | dualColumn[i]=tempModel[i]->dualColumnSolution(); |
2279 | | primalRow[i]=tempModel[i]->primalRowSolution(); |
2280 | | memcpy(dualColumn[i],model2->objective(), |
2281 | | numberColumns*sizeof(double)); |
2282 | | memset(primalRow[i], 0, numberRows * sizeof(double)); |
2283 | | tempModel[i]->clpMatrix()->transposeTimes(-1.0, |
2284 | | dualRow[i], |
2285 | | dualColumn[i]); |
2286 | | tempModel[i]->clpMatrix()->times(1.0, |
2287 | | primalColumn[i], |
2288 | | primalRow[i]); |
2289 | | tempModel[i]->checkSolutionInternal(); |
2290 | | printf("model %d - dual inf %g primal inf %g\n", |
2291 | | i,tempModel[i]->sumDualInfeasibilities(), |
2292 | | tempModel[i]->sumPrimalInfeasibilities()); |
2293 | | } |
2294 | | printf("What now\n"); |
2295 | | } else { |
2296 | | doubleIdiot=false; |
2297 | | } |
2298 | | } |
2299 | | if (!doubleIdiot) |
2300 | | info.crash(nPasses, model2->messageHandler(), model2->messagesPointer(), |
2301 | | (objective_->type() <2)); |
| 2247 | // move duals and just copy primal |
| 2248 | memcpy(temp.dualRowSolution(), dualModel2->primalColumnSolution(), |
| 2249 | numberRows * sizeof(double)); |
| 2250 | memcpy(temp.primalColumnSolution(), model2->primalColumnSolution(), |
| 2251 | numberColumns * sizeof(double)); |
| 2252 | #endif |
| 2253 | delete dualModel2; |
| 2254 | int numberRows = model2->numberRows(); |
| 2255 | int numberColumns = model2->numberColumns(); |
| 2256 | ClpSimplex *tempModel[2]; |
| 2257 | tempModel[0] = model2; |
| 2258 | tempModel[1] = &temp; |
| 2259 | const double *primalColumn[2]; |
| 2260 | const double *dualRow[2]; |
| 2261 | double *dualColumn[2]; |
| 2262 | double *primalRow[2]; |
| 2263 | for (int i = 0; i < 2; i++) { |
| 2264 | primalColumn[i] = tempModel[i]->primalColumnSolution(); |
| 2265 | dualRow[i] = tempModel[i]->dualRowSolution(); |
| 2266 | dualColumn[i] = tempModel[i]->dualColumnSolution(); |
| 2267 | primalRow[i] = tempModel[i]->primalRowSolution(); |
| 2268 | memcpy(dualColumn[i], model2->objective(), |
| 2269 | numberColumns * sizeof(double)); |
| 2270 | memset(primalRow[i], 0, numberRows * sizeof(double)); |
| 2271 | tempModel[i]->clpMatrix()->transposeTimes(-1.0, |
| 2272 | dualRow[i], |
| 2273 | dualColumn[i]); |
| 2274 | tempModel[i]->clpMatrix()->times(1.0, |
| 2275 | primalColumn[i], |
| 2276 | primalRow[i]); |
| 2277 | tempModel[i]->checkSolutionInternal(); |
| 2278 | printf("model %d - dual inf %g primal inf %g\n", |
| 2279 | i, tempModel[i]->sumDualInfeasibilities(), |
| 2280 | tempModel[i]->sumPrimalInfeasibilities()); |
| 2281 | } |
| 2282 | printf("What now\n"); |
| 2283 | } else { |
| 2284 | doubleIdiot = false; |
| 2285 | } |
| 2286 | } |
| 2287 | if (!doubleIdiot) |
| 2288 | info.crash(nPasses, model2->messageHandler(), model2->messagesPointer(), |
| 2289 | (objective_->type() < 2)); |
2303 | | info.crash(nPasses, model2->messageHandler(), model2->messagesPointer(),(objective_->type() <2)); |
2304 | | #endif |
2305 | | model2->scaling(saveScalingFlag); |
2306 | | #endif |
2307 | | time2 = CoinCpuTime(); |
2308 | | timeIdiot = time2 - timeX; |
2309 | | handler_->message(CLP_INTERVAL_TIMING, messages_) |
2310 | | << "Idiot Crash" << timeIdiot << time2 - time1 |
2311 | | << CoinMessageEol; |
2312 | | timeX = time2; |
2313 | | if (nPasses>100000&&nPasses<100500) { |
2314 | | // make sure no status left |
2315 | | model2->createStatus(); |
2316 | | // solve |
2317 | | if (model2->factorizationFrequency() == 200) { |
2318 | | // User did not touch preset |
2319 | | model2->defaultFactorizationFrequency(); |
2320 | | } |
2321 | | //int numberRows = model2->numberRows(); |
2322 | | int numberColumns = model2->numberColumns(); |
2323 | | // save duals |
2324 | | //double * saveDuals = CoinCopyOfArray(model2->dualRowSolution(),numberRows); |
2325 | | // for moment this only works on nug etc (i.e. all ==) |
2326 | | // needs row objective |
2327 | | double * saveObj = CoinCopyOfArray(model2->objective(),numberColumns); |
2328 | | double * pi = model2->dualRowSolution(); |
2329 | | model2->clpMatrix()->transposeTimes(-1.0, pi, model2->objective()); |
2330 | | // just primal values pass |
2331 | | double saveScale = model2->objectiveScale(); |
2332 | | model2->setObjectiveScale(1.0e-3); |
2333 | | model2->primal(2); |
2334 | | model2->writeMps("xx.mps"); |
2335 | | double * solution = model2->primalColumnSolution(); |
2336 | | double * upper = model2->columnUpper(); |
2337 | | for (int i=0;i<numberColumns;i++) { |
2338 | | if (solution[i]<100.0) |
2339 | | upper[i]=1000.0; |
2340 | | } |
2341 | | model2->setProblemStatus(-1); |
2342 | | model2->setObjectiveScale(saveScale); |
| 2291 | info.crash(nPasses, model2->messageHandler(), model2->messagesPointer(), (objective_->type() < 2)); |
| 2292 | #endif |
| 2293 | model2->scaling(saveScalingFlag); |
| 2294 | #endif |
| 2295 | time2 = CoinCpuTime(); |
| 2296 | timeIdiot = time2 - timeX; |
| 2297 | handler_->message(CLP_INTERVAL_TIMING, messages_) |
| 2298 | << "Idiot Crash" << timeIdiot << time2 - time1 |
| 2299 | << CoinMessageEol; |
| 2300 | timeX = time2; |
| 2301 | if (nPasses > 100000 && nPasses < 100500) { |
| 2302 | // make sure no status left |
| 2303 | model2->createStatus(); |
| 2304 | // solve |
| 2305 | if (model2->factorizationFrequency() == 200) { |
| 2306 | // User did not touch preset |
| 2307 | model2->defaultFactorizationFrequency(); |
| 2308 | } |
| 2309 | //int numberRows = model2->numberRows(); |
| 2310 | int numberColumns = model2->numberColumns(); |
| 2311 | // save duals |
| 2312 | //double * saveDuals = CoinCopyOfArray(model2->dualRowSolution(),numberRows); |
| 2313 | // for moment this only works on nug etc (i.e. all ==) |
| 2314 | // needs row objective |
| 2315 | double *saveObj = CoinCopyOfArray(model2->objective(), numberColumns); |
| 2316 | double *pi = model2->dualRowSolution(); |
| 2317 | model2->clpMatrix()->transposeTimes(-1.0, pi, model2->objective()); |
| 2318 | // just primal values pass |
| 2319 | double saveScale = model2->objectiveScale(); |
| 2320 | model2->setObjectiveScale(1.0e-3); |
| 2321 | model2->primal(2); |
| 2322 | model2->writeMps("xx.mps"); |
| 2323 | double *solution = model2->primalColumnSolution(); |
| 2324 | double *upper = model2->columnUpper(); |
| 2325 | for (int i = 0; i < numberColumns; i++) { |
| 2326 | if (solution[i] < 100.0) |
| 2327 | upper[i] = 1000.0; |
| 2328 | } |
| 2329 | model2->setProblemStatus(-1); |
| 2330 | model2->setObjectiveScale(saveScale); |
2469 | | // See if we have costed slacks |
2470 | | int * negSlack = new int[numberRows+numberColumns]; |
2471 | | int * posSlack = new int[numberRows+numberColumns]; |
2472 | | int * nextNegSlack = negSlack+numberRows; |
2473 | | int * nextPosSlack = posSlack+numberRows; |
2474 | | int iRow; |
2475 | | for (iRow = 0; iRow < numberRows+numberColumns; iRow++) { |
2476 | | negSlack[iRow] = -1; |
2477 | | posSlack[iRow] = -1; |
| 2456 | // See if we have costed slacks |
| 2457 | int *negSlack = new int[numberRows + numberColumns]; |
| 2458 | int *posSlack = new int[numberRows + numberColumns]; |
| 2459 | int *nextNegSlack = negSlack + numberRows; |
| 2460 | int *nextPosSlack = posSlack + numberRows; |
| 2461 | int iRow; |
| 2462 | for (iRow = 0; iRow < numberRows + numberColumns; iRow++) { |
| 2463 | negSlack[iRow] = -1; |
| 2464 | posSlack[iRow] = -1; |
| 2465 | } |
| 2466 | const double *element = model2->matrix()->getElements(); |
| 2467 | const int *row = model2->matrix()->getIndices(); |
| 2468 | const CoinBigIndex *columnStart = model2->matrix()->getVectorStarts(); |
| 2469 | const int *columnLength = model2->matrix()->getVectorLengths(); |
| 2470 | //bool allSlack = (numberRowsBasic==numberRows); |
| 2471 | for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) { |
| 2472 | if (!columnSolution[iColumn] || fabs(columnSolution[iColumn]) > 1.0e20) { |
| 2473 | double value = 0.0; |
| 2474 | if (columnLower[iColumn] > 0.0) |
| 2475 | value = columnLower[iColumn]; |
| 2476 | else if (columnUpper[iColumn] < 0.0) |
| 2477 | value = columnUpper[iColumn]; |
| 2478 | columnSolution[iColumn] = value; |
| 2479 | } |
| 2480 | if (columnLength[iColumn] == 1) { |
| 2481 | int jRow = row[columnStart[iColumn]]; |
| 2482 | if (!columnLower[iColumn]) { |
| 2483 | if (element[columnStart[iColumn]] > 0.0) { |
| 2484 | if (posSlack[jRow] < 0) { |
| 2485 | posSlack[jRow] = iColumn; |
| 2486 | } else { |
| 2487 | int jColumn = posSlack[jRow]; |
| 2488 | while (nextPosSlack[jColumn] >= 0) |
| 2489 | jColumn = nextPosSlack[jColumn]; |
| 2490 | nextPosSlack[jColumn] = iColumn; |
| 2491 | } |
| 2492 | } else if (element[columnStart[iColumn]] < 0.0) { |
| 2493 | if (negSlack[jRow] < 0) { |
| 2494 | negSlack[jRow] = iColumn; |
| 2495 | } else { |
| 2496 | int jColumn = negSlack[jRow]; |
| 2497 | while (nextNegSlack[jColumn] >= 0) |
| 2498 | jColumn = nextNegSlack[jColumn]; |
| 2499 | nextNegSlack[jColumn] = iColumn; |
| 2500 | } |
2479 | | const double * element = model2->matrix()->getElements(); |
2480 | | const int * row = model2->matrix()->getIndices(); |
2481 | | const CoinBigIndex * columnStart = model2->matrix()->getVectorStarts(); |
2482 | | const int * columnLength = model2->matrix()->getVectorLengths(); |
2483 | | //bool allSlack = (numberRowsBasic==numberRows); |
2484 | | for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) { |
2485 | | if (!columnSolution[iColumn] || fabs(columnSolution[iColumn]) > 1.0e20) { |
2486 | | double value = 0.0; |
2487 | | if (columnLower[iColumn] > 0.0) |
2488 | | value = columnLower[iColumn]; |
2489 | | else if (columnUpper[iColumn] < 0.0) |
2490 | | value = columnUpper[iColumn]; |
2491 | | columnSolution[iColumn] = value; |
2492 | | } |
2493 | | if (columnLength[iColumn] == 1) { |
2494 | | int jRow = row[columnStart[iColumn]]; |
2495 | | if (!columnLower[iColumn]) { |
2496 | | if (element[columnStart[iColumn]] > 0.0) { |
2497 | | if(posSlack[jRow] < 0) { |
2498 | | posSlack[jRow] = iColumn; |
2499 | | } else { |
2500 | | int jColumn=posSlack[jRow]; |
2501 | | while (nextPosSlack[jColumn]>=0) |
2502 | | jColumn=nextPosSlack[jColumn]; |
2503 | | nextPosSlack[jColumn]=iColumn; |
2504 | | } |
2505 | | } else if (element[columnStart[iColumn]] < 0.0) { |
2506 | | if(negSlack[jRow] < 0) { |
2507 | | negSlack[jRow] = iColumn; |
2508 | | } else { |
2509 | | int jColumn=negSlack[jRow]; |
2510 | | while (nextNegSlack[jColumn]>=0) |
2511 | | jColumn=nextNegSlack[jColumn]; |
2512 | | nextNegSlack[jColumn]=iColumn; |
2513 | | } |
2514 | | } |
2515 | | } else if (!columnUpper[iColumn]&& false) {// out for testing |
2516 | | if (element[columnStart[iColumn]] < 0.0 && posSlack[jRow] < 0) |
2517 | | posSlack[jRow] = iColumn; |
2518 | | else if (element[columnStart[iColumn]] > 0.0 && negSlack[jRow] < 0) |
2519 | | negSlack[jRow] = iColumn; |
2520 | | } |
2521 | | } |
| 2502 | } else if (!columnUpper[iColumn] && false) { // out for testing |
| 2503 | if (element[columnStart[iColumn]] < 0.0 && posSlack[jRow] < 0) |
| 2504 | posSlack[jRow] = iColumn; |
| 2505 | else if (element[columnStart[iColumn]] > 0.0 && negSlack[jRow] < 0) |
| 2506 | negSlack[jRow] = iColumn; |
| 2507 | } |
| 2508 | } |
| 2509 | } |
| 2510 | // now see what that does to row solution |
| 2511 | const double *objective = model2->objective(); |
| 2512 | double *rowSolution = model2->primalRowSolution(); |
| 2513 | CoinZeroN(rowSolution, numberRows); |
| 2514 | model2->clpMatrix()->times(1.0, columnSolution, rowSolution); |
| 2515 | // See if we can adjust using costed slacks |
| 2516 | double penalty = CoinMax(1.0e5, CoinMin(infeasibilityCost_ * 0.01, 1.0e10)) * optimizationDirection_; |
| 2517 | const double *lower = model2->rowLower(); |
| 2518 | const double *upper = model2->rowUpper(); |
| 2519 | for (iRow = 0; iRow < numberRows; iRow++) { |
| 2520 | if (lower[iRow] > rowSolution[iRow] + 1.0e-8) { |
| 2521 | int jColumn = posSlack[iRow]; |
| 2522 | if (jColumn >= 0) { |
| 2523 | // sort if more than one |
| 2524 | int nPos = 1; |
| 2525 | sort[0] = jColumn; |
| 2526 | weight[0] = objective[jColumn]; |
| 2527 | while (nextPosSlack[jColumn] >= 0) { |
| 2528 | jColumn = nextPosSlack[jColumn]; |
| 2529 | sort[nPos] = jColumn; |
| 2530 | weight[nPos++] = objective[jColumn]; |
2523 | | // now see what that does to row solution |
2524 | | const double * objective = model2->objective(); |
2525 | | double * rowSolution = model2->primalRowSolution(); |
2526 | | CoinZeroN (rowSolution, numberRows); |
2527 | | model2->clpMatrix()->times(1.0, columnSolution, rowSolution); |
2528 | | // See if we can adjust using costed slacks |
2529 | | double penalty = CoinMax(1.0e5, CoinMin(infeasibilityCost_ * 0.01, 1.0e10)) * optimizationDirection_; |
2530 | | const double * lower = model2->rowLower(); |
2531 | | const double * upper = model2->rowUpper(); |
2532 | | for (iRow = 0; iRow < numberRows; iRow++) { |
2533 | | if (lower[iRow] > rowSolution[iRow] + 1.0e-8) { |
2534 | | int jColumn = posSlack[iRow]; |
2535 | | if (jColumn >= 0) { |
2536 | | // sort if more than one |
2537 | | int nPos=1; |
2538 | | sort[0]=jColumn; |
2539 | | weight[0]=objective[jColumn]; |
2540 | | while (nextPosSlack[jColumn]>=0) { |
2541 | | jColumn=nextPosSlack[jColumn]; |
2542 | | sort[nPos]=jColumn; |
2543 | | weight[nPos++]=objective[jColumn]; |
2544 | | } |
2545 | | if (nPos>1) { |
2546 | | CoinSort_2(weight,weight+nPos,sort); |
2547 | | for (int i=0;i<nPos;i++) { |
2548 | | double difference = lower[iRow] - rowSolution[iRow]; |
2549 | | jColumn=sort[i]; |
2550 | | double elementValue = element[columnStart[jColumn]]; |
2551 | | assert (elementValue > 0.0); |
2552 | | double value=columnSolution[jColumn]; |
2553 | | double movement = CoinMin(difference / elementValue, columnUpper[jColumn]-value); |
2554 | | columnSolution[jColumn] += movement; |
2555 | | rowSolution[iRow] += movement * elementValue; |
2556 | | } |
2557 | | continue; |
2558 | | } |
2559 | | if (jColumn<0||columnSolution[jColumn]) |
2560 | | continue; |
2561 | | double difference = lower[iRow] - rowSolution[iRow]; |
2562 | | double elementValue = element[columnStart[jColumn]]; |
2563 | | if (elementValue > 0.0) { |
2564 | | double movement = CoinMin(difference / elementValue, columnUpper[jColumn]); |
2565 | | columnSolution[jColumn] = movement; |
2566 | | rowSolution[iRow] += movement * elementValue; |
2567 | | } else { |
2568 | | double movement = CoinMax(difference / elementValue, columnLower[jColumn]); |
2569 | | columnSolution[jColumn] = movement; |
2570 | | rowSolution[iRow] += movement * elementValue; |
2571 | | } |
2572 | | } |
2573 | | } else if (upper[iRow] < rowSolution[iRow] - 1.0e-8) { |
2574 | | int jColumn = negSlack[iRow]; |
2575 | | if (jColumn >= 0) { |
2576 | | // sort if more than one |
2577 | | int nNeg=1; |
2578 | | sort[0]=jColumn; |
2579 | | weight[0]=objective[jColumn]; |
2580 | | while (nextNegSlack[jColumn]>=0) { |
2581 | | jColumn=nextNegSlack[jColumn]; |
2582 | | sort[nNeg]=jColumn; |
2583 | | weight[nNeg++]=objective[jColumn]; |
2584 | | } |
2585 | | if (nNeg>1) { |
2586 | | CoinSort_2(weight,weight+nNeg,sort); |
2587 | | for (int i=0;i<nNeg;i++) { |
2588 | | double difference = rowSolution[iRow]-upper[iRow]; |
2589 | | jColumn=sort[i]; |
2590 | | double elementValue = element[columnStart[jColumn]]; |
2591 | | assert (elementValue < 0.0); |
2592 | | double value=columnSolution[jColumn]; |
2593 | | double movement = CoinMin(difference / -elementValue, columnUpper[jColumn]-value); |
2594 | | columnSolution[jColumn] += movement; |
2595 | | rowSolution[iRow] += movement * elementValue; |
2596 | | } |
2597 | | continue; |
2598 | | } |
2599 | | if (jColumn<0||columnSolution[jColumn]) |
2600 | | continue; |
2601 | | double difference = upper[iRow] - rowSolution[iRow]; |
2602 | | double elementValue = element[columnStart[jColumn]]; |
2603 | | if (elementValue < 0.0) { |
2604 | | double movement = CoinMin(difference / elementValue, columnUpper[jColumn]); |
2605 | | columnSolution[jColumn] = movement; |
2606 | | rowSolution[iRow] += movement * elementValue; |
2607 | | } else { |
2608 | | double movement = CoinMax(difference / elementValue, columnLower[jColumn]); |
2609 | | columnSolution[jColumn] = movement; |
2610 | | rowSolution[iRow] += movement * elementValue; |
2611 | | } |
2612 | | } |
2613 | | } |
| 2532 | if (nPos > 1) { |
| 2533 | CoinSort_2(weight, weight + nPos, sort); |
| 2534 | for (int i = 0; i < nPos; i++) { |
| 2535 | double difference = lower[iRow] - rowSolution[iRow]; |
| 2536 | jColumn = sort[i]; |
| 2537 | double elementValue = element[columnStart[jColumn]]; |
| 2538 | assert(elementValue > 0.0); |
| 2539 | double value = columnSolution[jColumn]; |
| 2540 | double movement = CoinMin(difference / elementValue, columnUpper[jColumn] - value); |
| 2541 | columnSolution[jColumn] += movement; |
| 2542 | rowSolution[iRow] += movement * elementValue; |
| 2543 | } |
| 2544 | continue; |
2684 | | double * saveLower = NULL; |
2685 | | double * saveUpper = NULL; |
2686 | | if (largest < -2.01 * smallest) { // switch off for now |
2687 | | // perturb - so switch off standard |
2688 | | model2->setPerturbation(100); |
2689 | | saveLower = new double[numberRows]; |
2690 | | CoinMemcpyN(model2->rowLower_, numberRows, saveLower); |
2691 | | saveUpper = new double[numberRows]; |
2692 | | CoinMemcpyN(model2->rowUpper_, numberRows, saveUpper); |
2693 | | double * lower = model2->rowLower(); |
2694 | | double * upper = model2->rowUpper(); |
2695 | | for (iRow = 0; iRow < numberRows; iRow++) { |
2696 | | double lowerValue = lower[iRow], upperValue = upper[iRow]; |
2697 | | double value = randomNumberGenerator_.randomDouble(); |
2698 | | if (upperValue > lowerValue + primalTolerance_) { |
2699 | | if (lowerValue > -1.0e20 && lowerValue) |
2700 | | lowerValue -= value * 1.0e-4 * fabs(lowerValue); |
2701 | | if (upperValue < 1.0e20 && upperValue) |
2702 | | upperValue += value * 1.0e-4 * fabs(upperValue); |
2703 | | } else if (upperValue > 0.0) { |
2704 | | upperValue -= value * 1.0e-4 * fabs(lowerValue); |
2705 | | lowerValue -= value * 1.0e-4 * fabs(lowerValue); |
2706 | | } else if (upperValue < 0.0) { |
2707 | | upperValue += value * 1.0e-4 * fabs(lowerValue); |
2708 | | lowerValue += value * 1.0e-4 * fabs(lowerValue); |
2709 | | } else { |
2710 | | } |
2711 | | lower[iRow] = lowerValue; |
2712 | | upper[iRow] = upperValue; |
2713 | | } |
2714 | | } |
2715 | | int i; |
2716 | | // Just do this number of passes in Sprint |
2717 | | if (doSprint > 0) |
2718 | | maxSprintPass = options.getExtraInfo(1); |
2719 | | // but if big use to get ratio |
2720 | | double ratio = 3; |
| 2599 | } |
| 2600 | } |
| 2601 | } |
| 2602 | delete[] negSlack; |
| 2603 | delete[] posSlack; |
| 2604 | int nRow = numberRows; |
| 2605 | bool network = false; |
| 2606 | if (dynamic_cast< ClpNetworkMatrix * >(matrix_)) { |
| 2607 | network = true; |
| 2608 | nRow *= 2; |
| 2609 | } |
| 2610 | CoinBigIndex *addStarts = new CoinBigIndex[nRow + 1]; |
| 2611 | int *addRow = new int[nRow]; |
| 2612 | double *addElement = new double[nRow]; |
| 2613 | addStarts[0] = 0; |
| 2614 | int numberArtificials = 0; |
| 2615 | int numberAdd = 0; |
| 2616 | double *addCost = new double[numberRows]; |
| 2617 | for (iRow = 0; iRow < numberRows; iRow++) { |
| 2618 | if (lower[iRow] > rowSolution[iRow] + 1.0e-8) { |
| 2619 | addRow[numberAdd] = iRow; |
| 2620 | addElement[numberAdd++] = 1.0; |
| 2621 | if (network) { |
| 2622 | addRow[numberAdd] = numberRows; |
| 2623 | addElement[numberAdd++] = -1.0; |
| 2624 | } |
| 2625 | addCost[numberArtificials] = penalty; |
| 2626 | numberArtificials++; |
| 2627 | addStarts[numberArtificials] = numberAdd; |
| 2628 | } else if (upper[iRow] < rowSolution[iRow] - 1.0e-8) { |
| 2629 | addRow[numberAdd] = iRow; |
| 2630 | addElement[numberAdd++] = -1.0; |
| 2631 | if (network) { |
| 2632 | addRow[numberAdd] = numberRows; |
| 2633 | addElement[numberAdd++] = 1.0; |
| 2634 | } |
| 2635 | addCost[numberArtificials] = penalty; |
| 2636 | numberArtificials++; |
| 2637 | addStarts[numberArtificials] = numberAdd; |
| 2638 | } |
| 2639 | } |
| 2640 | if (numberArtificials) { |
| 2641 | // need copy so as not to disturb original |
| 2642 | model2 = new ClpSimplex(*model2); |
| 2643 | if (network) { |
| 2644 | // network - add a null row |
| 2645 | model2->addRow(0, NULL, NULL, -COIN_DBL_MAX, COIN_DBL_MAX); |
| 2646 | numberRows++; |
| 2647 | } |
| 2648 | model2->addColumns(numberArtificials, NULL, NULL, addCost, |
| 2649 | addStarts, addRow, addElement); |
| 2650 | } |
| 2651 | delete[] addStarts; |
| 2652 | delete[] addRow; |
| 2653 | delete[] addElement; |
| 2654 | delete[] addCost; |
| 2655 | // look at rhs to see if to perturb |
| 2656 | double largest = 0.0; |
| 2657 | double smallest = 1.0e30; |
| 2658 | for (iRow = 0; iRow < numberRows; iRow++) { |
| 2659 | double value; |
| 2660 | value = fabs(model2->rowLower_[iRow]); |
| 2661 | if (value && value < 1.0e30) { |
| 2662 | largest = CoinMax(largest, value); |
| 2663 | smallest = CoinMin(smallest, value); |
| 2664 | } |
| 2665 | value = fabs(model2->rowUpper_[iRow]); |
| 2666 | if (value && value < 1.0e30) { |
| 2667 | largest = CoinMax(largest, value); |
| 2668 | smallest = CoinMin(smallest, value); |
| 2669 | } |
| 2670 | } |
| 2671 | double *saveLower = NULL; |
| 2672 | double *saveUpper = NULL; |
| 2673 | if (largest < -2.01 * smallest) { // switch off for now |
| 2674 | // perturb - so switch off standard |
| 2675 | model2->setPerturbation(100); |
| 2676 | saveLower = new double[numberRows]; |
| 2677 | CoinMemcpyN(model2->rowLower_, numberRows, saveLower); |
| 2678 | saveUpper = new double[numberRows]; |
| 2679 | CoinMemcpyN(model2->rowUpper_, numberRows, saveUpper); |
| 2680 | double *lower = model2->rowLower(); |
| 2681 | double *upper = model2->rowUpper(); |
| 2682 | for (iRow = 0; iRow < numberRows; iRow++) { |
| 2683 | double lowerValue = lower[iRow], upperValue = upper[iRow]; |
| 2684 | double value = randomNumberGenerator_.randomDouble(); |
| 2685 | if (upperValue > lowerValue + primalTolerance_) { |
| 2686 | if (lowerValue > -1.0e20 && lowerValue) |
| 2687 | lowerValue -= value * 1.0e-4 * fabs(lowerValue); |
| 2688 | if (upperValue < 1.0e20 && upperValue) |
| 2689 | upperValue += value * 1.0e-4 * fabs(upperValue); |
| 2690 | } else if (upperValue > 0.0) { |
| 2691 | upperValue -= value * 1.0e-4 * fabs(lowerValue); |
| 2692 | lowerValue -= value * 1.0e-4 * fabs(lowerValue); |
| 2693 | } else if (upperValue < 0.0) { |
| 2694 | upperValue += value * 1.0e-4 * fabs(lowerValue); |
| 2695 | lowerValue += value * 1.0e-4 * fabs(lowerValue); |
| 2696 | } else { |
| 2697 | } |
| 2698 | lower[iRow] = lowerValue; |
| 2699 | upper[iRow] = upperValue; |
| 2700 | } |
| 2701 | } |
| 2702 | int i; |
| 2703 | // Just do this number of passes in Sprint |
| 2704 | if (doSprint > 0) |
| 2705 | maxSprintPass = options.getExtraInfo(1); |
| 2706 | // but if big use to get ratio |
| 2707 | double ratio = 3; |
2730 | | printf("%d passes wanted with ratio of %g\n", maxSprintPass, ratio); |
2731 | | #endif |
2732 | | } |
2733 | | // Just take this number of columns in small problem |
2734 | | int smallNumberColumns = static_cast<int> (CoinMin(ratio * numberRows, static_cast<double> (numberColumns))); |
2735 | | smallNumberColumns = CoinMax(smallNumberColumns, 3000); |
2736 | | smallNumberColumns = CoinMin(smallNumberColumns, numberColumns); |
2737 | | int saveSmallNumber = smallNumberColumns; |
2738 | | bool emergencyMode=false; |
2739 | | //int smallNumberColumns = CoinMin(12*numberRows/10,numberColumns); |
2740 | | //smallNumberColumns = CoinMax(smallNumberColumns,3000); |
2741 | | //smallNumberColumns = CoinMax(smallNumberColumns,numberRows+1000); |
2742 | | // redo as may have changed |
2743 | | columnLower = model2->columnLower(); |
2744 | | columnUpper = model2->columnUpper(); |
2745 | | columnSolution = model2->primalColumnSolution(); |
2746 | | // Set up initial list |
2747 | | numberSort = 0; |
2748 | | if (numberArtificials) { |
2749 | | numberSort = numberArtificials; |
2750 | | for (i = 0; i < numberSort; i++) |
2751 | | sort[i] = i + originalNumberColumns; |
2752 | | } |
2753 | | // put in free |
2754 | | // maybe a solution there already |
2755 | | for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) { |
2756 | | if (model2->getColumnStatus(iColumn) == basic|| |
2757 | | (columnLower[iColumn]<-1.0e30&&columnUpper[iColumn]>1.0e30)) |
2758 | | sort[numberSort++] = iColumn; |
2759 | | } |
2760 | | for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) { |
2761 | | if (model2->getColumnStatus(iColumn) != basic) { |
2762 | | if (columnSolution[iColumn] > columnLower[iColumn] && |
2763 | | columnSolution[iColumn] < columnUpper[iColumn] && |
2764 | | columnSolution[iColumn]) |
2765 | | sort[numberSort++] = iColumn; |
2766 | | } |
2767 | | } |
2768 | | numberSort = CoinMin(numberSort, smallNumberColumns); |
| 2717 | printf("%d passes wanted with ratio of %g\n", maxSprintPass, ratio); |
| 2718 | #endif |
| 2719 | } |
| 2720 | // Just take this number of columns in small problem |
| 2721 | int smallNumberColumns = static_cast< int >(CoinMin(ratio * numberRows, static_cast< double >(numberColumns))); |
| 2722 | smallNumberColumns = CoinMax(smallNumberColumns, 3000); |
| 2723 | smallNumberColumns = CoinMin(smallNumberColumns, numberColumns); |
| 2724 | int saveSmallNumber = smallNumberColumns; |
| 2725 | bool emergencyMode = false; |
| 2726 | //int smallNumberColumns = CoinMin(12*numberRows/10,numberColumns); |
| 2727 | //smallNumberColumns = CoinMax(smallNumberColumns,3000); |
| 2728 | //smallNumberColumns = CoinMax(smallNumberColumns,numberRows+1000); |
| 2729 | // redo as may have changed |
| 2730 | columnLower = model2->columnLower(); |
| 2731 | columnUpper = model2->columnUpper(); |
| 2732 | columnSolution = model2->primalColumnSolution(); |
| 2733 | // Set up initial list |
| 2734 | numberSort = 0; |
| 2735 | if (numberArtificials) { |
| 2736 | numberSort = numberArtificials; |
| 2737 | for (i = 0; i < numberSort; i++) |
| 2738 | sort[i] = i + originalNumberColumns; |
| 2739 | } |
| 2740 | // put in free |
| 2741 | // maybe a solution there already |
| 2742 | for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) { |
| 2743 | if (model2->getColumnStatus(iColumn) == basic || (columnLower[iColumn] < -1.0e30 && columnUpper[iColumn] > 1.0e30)) |
| 2744 | sort[numberSort++] = iColumn; |
| 2745 | } |
| 2746 | for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) { |
| 2747 | if (model2->getColumnStatus(iColumn) != basic) { |
| 2748 | if (columnSolution[iColumn] > columnLower[iColumn] && columnSolution[iColumn] < columnUpper[iColumn] && columnSolution[iColumn]) |
| 2749 | sort[numberSort++] = iColumn; |
| 2750 | } |
| 2751 | } |
| 2752 | numberSort = CoinMin(numberSort, smallNumberColumns); |
2774 | | int iPass; |
2775 | | double lastObjective[] = {1.0e31,1.0e31}; |
2776 | | // It will be safe to allow dense |
2777 | | model2->setInitialDenseFactorization(true); |
2778 | | |
2779 | | // We will be using all rows |
2780 | | int * whichRows = new int [numberRows]; |
2781 | | for (iRow = 0; iRow < numberRows; iRow++) |
2782 | | whichRows[iRow] = iRow; |
2783 | | double originalOffset; |
2784 | | model2->getDblParam(ClpObjOffset, originalOffset); |
2785 | | int totalIterations = 0; |
2786 | | double lastSumArtificials = COIN_DBL_MAX; |
2787 | | int originalMaxSprintPass = maxSprintPass; |
2788 | | maxSprintPass = 20; // so we do that many if infeasible |
2789 | | for (iPass = 0; iPass < maxSprintPass; iPass++) { |
2790 | | //printf("Bug until submodel new version\n"); |
2791 | | //CoinSort_2(sort,sort+numberSort,weight); |
2792 | | // Create small problem |
2793 | | ClpSimplex small(model2, numberRows, whichRows, numberSort, sort); |
2794 | | small.setPerturbation(model2->perturbation()); |
2795 | | small.setInfeasibilityCost(model2->infeasibilityCost()); |
2796 | | if (model2->factorizationFrequency() == 200) { |
2797 | | // User did not touch preset |
2798 | | small.defaultFactorizationFrequency(); |
2799 | | } |
2800 | | // now see what variables left out do to row solution |
2801 | | double * rowSolution = model2->primalRowSolution(); |
2802 | | double * sumFixed = new double[numberRows]; |
2803 | | CoinZeroN (sumFixed, numberRows); |
2804 | | int iRow, iColumn; |
2805 | | // zero out ones in small problem |
2806 | | for (iColumn = 0; iColumn < numberSort; iColumn++) { |
2807 | | int kColumn = sort[iColumn]; |
2808 | | fullSolution[kColumn] = 0.0; |
2809 | | } |
2810 | | // Get objective offset |
2811 | | const double * objective = model2->objective(); |
2812 | | double offset = 0.0; |
2813 | | for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) |
2814 | | offset += fullSolution[iColumn] * objective[iColumn]; |
| 2762 | // We will be using all rows |
| 2763 | int *whichRows = new int[numberRows]; |
| 2764 | for (iRow = 0; iRow < numberRows; iRow++) |
| 2765 | whichRows[iRow] = iRow; |
| 2766 | double originalOffset; |
| 2767 | model2->getDblParam(ClpObjOffset, originalOffset); |
| 2768 | int totalIterations = 0; |
| 2769 | double lastSumArtificials = COIN_DBL_MAX; |
| 2770 | int originalMaxSprintPass = maxSprintPass; |
| 2771 | maxSprintPass = 20; // so we do that many if infeasible |
| 2772 | for (iPass = 0; iPass < maxSprintPass; iPass++) { |
| 2773 | //printf("Bug until submodel new version\n"); |
| 2774 | //CoinSort_2(sort,sort+numberSort,weight); |
| 2775 | // Create small problem |
| 2776 | ClpSimplex small(model2, numberRows, whichRows, numberSort, sort); |
| 2777 | small.setPerturbation(model2->perturbation()); |
| 2778 | small.setInfeasibilityCost(model2->infeasibilityCost()); |
| 2779 | if (model2->factorizationFrequency() == 200) { |
| 2780 | // User did not touch preset |
| 2781 | small.defaultFactorizationFrequency(); |
| 2782 | } |
| 2783 | // now see what variables left out do to row solution |
| 2784 | double *rowSolution = model2->primalRowSolution(); |
| 2785 | double *sumFixed = new double[numberRows]; |
| 2786 | CoinZeroN(sumFixed, numberRows); |
| 2787 | int iRow, iColumn; |
| 2788 | // zero out ones in small problem |
| 2789 | for (iColumn = 0; iColumn < numberSort; iColumn++) { |
| 2790 | int kColumn = sort[iColumn]; |
| 2791 | fullSolution[kColumn] = 0.0; |
| 2792 | } |
| 2793 | // Get objective offset |
| 2794 | const double *objective = model2->objective(); |
| 2795 | double offset = 0.0; |
| 2796 | for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) |
| 2797 | offset += fullSolution[iColumn] * objective[iColumn]; |
2827 | | double * lower = small.rowLower(); |
2828 | | double * upper = small.rowUpper(); |
2829 | | for (iRow = 0; iRow < numberRows; iRow++) { |
2830 | | if (lower[iRow] > -1.0e50) |
2831 | | lower[iRow] -= sumFixed[iRow]; |
2832 | | if (upper[iRow] < 1.0e50) |
2833 | | upper[iRow] -= sumFixed[iRow]; |
2834 | | rowSolution[iRow] -= sumFixed[iRow]; |
2835 | | } |
2836 | | delete [] sumFixed; |
2837 | | // Solve |
2838 | | if (interrupt) |
2839 | | currentModel = &small; |
2840 | | small.defaultFactorizationFrequency(); |
2841 | | if (emergencyMode) { |
2842 | | // not much happening so big model |
2843 | | int options=small.moreSpecialOptions(); |
2844 | | small.setMoreSpecialOptions(options|1048576); |
2845 | | small.setMaximumIterations(1000400); |
2846 | | small.setPerturbation(100); |
2847 | | } |
2848 | | if (dynamic_cast< ClpPackedMatrix*>(matrix_)) { |
2849 | | // See if original wanted vector |
2850 | | ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_); |
2851 | | ClpMatrixBase * matrix = small.clpMatrix(); |
2852 | | if (dynamic_cast< ClpPackedMatrix*>(matrix) && clpMatrixO->wantsSpecialColumnCopy()) { |
2853 | | ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix); |
2854 | | clpMatrix->makeSpecialColumnCopy(); |
2855 | | small.primal(1); |
2856 | | clpMatrix->releaseSpecialColumnCopy(); |
2857 | | } else { |
| 2810 | double *lower = small.rowLower(); |
| 2811 | double *upper = small.rowUpper(); |
| 2812 | for (iRow = 0; iRow < numberRows; iRow++) { |
| 2813 | if (lower[iRow] > -1.0e50) |
| 2814 | lower[iRow] -= sumFixed[iRow]; |
| 2815 | if (upper[iRow] < 1.0e50) |
| 2816 | upper[iRow] -= sumFixed[iRow]; |
| 2817 | rowSolution[iRow] -= sumFixed[iRow]; |
| 2818 | } |
| 2819 | delete[] sumFixed; |
| 2820 | // Solve |
| 2821 | if (interrupt) |
| 2822 | currentModel = &small; |
| 2823 | small.defaultFactorizationFrequency(); |
| 2824 | if (emergencyMode) { |
| 2825 | // not much happening so big model |
| 2826 | int options = small.moreSpecialOptions(); |
| 2827 | small.setMoreSpecialOptions(options | 1048576); |
| 2828 | small.setMaximumIterations(1000400); |
| 2829 | small.setPerturbation(100); |
| 2830 | } |
| 2831 | if (dynamic_cast< ClpPackedMatrix * >(matrix_)) { |
| 2832 | // See if original wanted vector |
| 2833 | ClpPackedMatrix *clpMatrixO = dynamic_cast< ClpPackedMatrix * >(matrix_); |
| 2834 | ClpMatrixBase *matrix = small.clpMatrix(); |
| 2835 | if (dynamic_cast< ClpPackedMatrix * >(matrix) && clpMatrixO->wantsSpecialColumnCopy()) { |
| 2836 | ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(matrix); |
| 2837 | clpMatrix->makeSpecialColumnCopy(); |
| 2838 | small.primal(1); |
| 2839 | clpMatrix->releaseSpecialColumnCopy(); |
| 2840 | } else { |
2866 | | if (iPass||!numberArtificials) |
2867 | | small.primal(1); |
2868 | | else |
2869 | | small.dual(0); |
2870 | | #endif |
2871 | | if (emergencyMode) { |
2872 | | if (small.problemStatus()==3) |
2873 | | small.setProblemStatus(0); |
2874 | | smallNumberColumns=saveSmallNumber; |
2875 | | emergencyMode=false; |
2876 | | double * temp = new double [numberRows]; |
2877 | | memset(temp,0,numberRows*sizeof(double)); |
2878 | | double * solution = small.primalColumnSolution(); |
2879 | | small.matrix()->times(solution,temp); |
2880 | | double sumInf=0.0; |
2881 | | double * lower = small.rowLower(); |
2882 | | double * upper = small.rowUpper(); |
2883 | | for (int iRow=0;iRow<numberRows;iRow++) { |
2884 | | if (temp[iRow]>upper[iRow]) |
2885 | | sumInf += temp[iRow]-upper[iRow]; |
2886 | | else if (temp[iRow]<lower[iRow]) |
2887 | | sumInf += lower[iRow]-temp[iRow]; |
2888 | | } |
2889 | | printf("row inf %g\n",sumInf); |
2890 | | sumInf=0.0; |
2891 | | lower = small.columnLower(); |
2892 | | upper = small.columnUpper(); |
2893 | | for (int iColumn=0;iColumn<small.numberColumns();iColumn++) { |
2894 | | if (solution[iColumn]>upper[iColumn]) |
2895 | | sumInf += solution[iColumn]-upper[iColumn]; |
2896 | | else if (solution[iColumn]<lower[iColumn]) |
2897 | | sumInf += lower[iColumn]-solution[iColumn]; |
2898 | | } |
2899 | | printf("column inf %g\n",sumInf); |
2900 | | delete [] temp; |
2901 | | } |
2902 | | if (small.problemStatus()) { |
2903 | | int numberIterations=small.numberIterations(); |
2904 | | small.dual(0); |
2905 | | small.setNumberIterations(small.numberIterations()+numberIterations); |
2906 | | } |
| 2849 | if (iPass || !numberArtificials) |
| 2850 | small.primal(1); |
| 2851 | else |
| 2852 | small.dual(0); |
| 2853 | #endif |
| 2854 | if (emergencyMode) { |
| 2855 | if (small.problemStatus() == 3) |
| 2856 | small.setProblemStatus(0); |
| 2857 | smallNumberColumns = saveSmallNumber; |
| 2858 | emergencyMode = false; |
| 2859 | double *temp = new double[numberRows]; |
| 2860 | memset(temp, 0, numberRows * sizeof(double)); |
| 2861 | double *solution = small.primalColumnSolution(); |
| 2862 | small.matrix()->times(solution, temp); |
| 2863 | double sumInf = 0.0; |
| 2864 | double *lower = small.rowLower(); |
| 2865 | double *upper = small.rowUpper(); |
| 2866 | for (int iRow = 0; iRow < numberRows; iRow++) { |
| 2867 | if (temp[iRow] > upper[iRow]) |
| 2868 | sumInf += temp[iRow] - upper[iRow]; |
| 2869 | else if (temp[iRow] < lower[iRow]) |
| 2870 | sumInf += lower[iRow] - temp[iRow]; |
| 2871 | } |
| 2872 | printf("row inf %g\n", sumInf); |
| 2873 | sumInf = 0.0; |
| 2874 | lower = small.columnLower(); |
| 2875 | upper = small.columnUpper(); |
| 2876 | for (int iColumn = 0; iColumn < small.numberColumns(); iColumn++) { |
| 2877 | if (solution[iColumn] > upper[iColumn]) |
| 2878 | sumInf += solution[iColumn] - upper[iColumn]; |
| 2879 | else if (solution[iColumn] < lower[iColumn]) |
| 2880 | sumInf += lower[iColumn] - solution[iColumn]; |
| 2881 | } |
| 2882 | printf("column inf %g\n", sumInf); |
| 2883 | delete[] temp; |
| 2884 | } |
| 2885 | if (small.problemStatus()) { |
| 2886 | int numberIterations = small.numberIterations(); |
| 2887 | small.dual(0); |
| 2888 | small.setNumberIterations(small.numberIterations() + numberIterations); |
| 2889 | } |
2930 | | small.primal(1); |
2931 | | #endif |
2932 | | if (small2->problemStatus()) |
2933 | | small.primal(1); |
2934 | | } |
2935 | | delete small2; |
2936 | | } else { |
2937 | | small.primal(1); |
2938 | | } |
2939 | | delete [] whichRow; |
2940 | | delete [] whichColumn; |
2941 | | #endif |
2942 | | } |
2943 | | } else { |
2944 | | small.primal(1); |
2945 | | } |
2946 | | int smallIterations=small.numberIterations(); |
2947 | | totalIterations += smallIterations; |
2948 | | if (2*smallIterations<CoinMin(numberRows,1000) && iPass) { |
2949 | | int oldNumber=smallNumberColumns; |
2950 | | if (smallIterations<100) |
2951 | | smallNumberColumns *= 1.2; |
2952 | | else |
2953 | | smallNumberColumns *= 1.1; |
2954 | | smallNumberColumns = CoinMin(smallNumberColumns, numberColumns); |
2955 | | if (smallIterations<200 && smallNumberColumns>2*saveSmallNumber) { |
2956 | | // try kicking it |
2957 | | emergencyMode=true; |
2958 | | smallNumberColumns=numberColumns; |
2959 | | } |
2960 | | // smallNumberColumns = CoinMin(smallNumberColumns, 3*saveSmallNumber); |
2961 | | char line[100]; |
2962 | | sprintf(line,"sample size increased from %d to %d", |
2963 | | oldNumber,smallNumberColumns); |
2964 | | handler_->message(CLP_GENERAL, messages_) |
2965 | | << line |
2966 | | << CoinMessageEol; |
2967 | | } |
2968 | | // move solution back |
2969 | | const double * solution = small.primalColumnSolution(); |
2970 | | for (iColumn = 0; iColumn < numberSort; iColumn++) { |
2971 | | int kColumn = sort[iColumn]; |
2972 | | model2->setColumnStatus(kColumn, small.getColumnStatus(iColumn)); |
2973 | | fullSolution[kColumn] = solution[iColumn]; |
2974 | | } |
2975 | | for (iRow = 0; iRow < numberRows; iRow++) |
2976 | | model2->setRowStatus(iRow, small.getRowStatus(iRow)); |
2977 | | CoinMemcpyN(small.primalRowSolution(), |
2978 | | numberRows, model2->primalRowSolution()); |
2979 | | double sumArtificials = 0.0; |
2980 | | for (i = 0; i < numberArtificials; i++) |
2981 | | sumArtificials += fullSolution[i + originalNumberColumns]; |
2982 | | if (sumArtificials && iPass > 5 && sumArtificials >= lastSumArtificials) { |
2983 | | // increase costs |
2984 | | double * cost = model2->objective() + originalNumberColumns; |
2985 | | double newCost = CoinMin(1.0e10, cost[0] * 1.5); |
2986 | | for (i = 0; i < numberArtificials; i++) |
2987 | | cost[i] = newCost; |
2988 | | } |
2989 | | lastSumArtificials = sumArtificials; |
2990 | | // get reduced cost for large problem |
2991 | | double * djs = model2->dualColumnSolution(); |
2992 | | CoinMemcpyN(model2->objective(), numberColumns, djs); |
2993 | | model2->clpMatrix()->transposeTimes(-1.0, small.dualRowSolution(), djs); |
2994 | | int numberNegative = 0; |
2995 | | double sumNegative = 0.0; |
2996 | | // now massage weight so all basic in plus good djs |
2997 | | // first count and do basic |
2998 | | numberSort = 0; |
2999 | | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
3000 | | double dj = djs[iColumn] * optimizationDirection_; |
3001 | | double value = fullSolution[iColumn]; |
3002 | | if (model2->getColumnStatus(iColumn) == ClpSimplex::basic) { |
3003 | | sort[numberSort++] = iColumn; |
3004 | | } else if (dj < -dualTolerance_ && value < columnUpper[iColumn]) { |
3005 | | numberNegative++; |
3006 | | sumNegative -= dj; |
3007 | | } else if (dj > dualTolerance_ && value > columnLower[iColumn]) { |
3008 | | numberNegative++; |
3009 | | sumNegative += dj; |
3010 | | } |
3011 | | } |
3012 | | handler_->message(CLP_SPRINT, messages_) |
3013 | | << iPass + 1 << small.numberIterations() << small.objectiveValue() << sumNegative |
3014 | | << numberNegative |
3015 | | << CoinMessageEol; |
3016 | | if (sumArtificials < 1.0e-8 && originalMaxSprintPass >= 0) { |
3017 | | maxSprintPass = iPass + originalMaxSprintPass; |
3018 | | originalMaxSprintPass = -1; |
3019 | | } |
3020 | | if (iPass > 20) |
3021 | | sumArtificials = 0.0; |
3022 | | if ((small.objectiveValue()*optimizationDirection_ > lastObjective[1] - 1.0e-7 && iPass > 15 && sumArtificials < 1.0e-8 && maxSprintPass<200) || |
3023 | | (!small.numberIterations() && iPass) || |
3024 | | iPass == maxSprintPass - 1 || small.status() == 3) { |
| 2912 | small.primal(1); |
| 2913 | #endif |
| 2914 | if (small2->problemStatus()) |
| 2915 | small.primal(1); |
| 2916 | } |
| 2917 | delete small2; |
| 2918 | } else { |
| 2919 | small.primal(1); |
| 2920 | } |
| 2921 | delete[] whichRow; |
| 2922 | delete[] whichColumn; |
| 2923 | #endif |
| 2924 | } |
| 2925 | } else { |
| 2926 | small.primal(1); |
| 2927 | } |
| 2928 | int smallIterations = small.numberIterations(); |
| 2929 | totalIterations += smallIterations; |
| 2930 | if (2 * smallIterations < CoinMin(numberRows, 1000) && iPass) { |
| 2931 | int oldNumber = smallNumberColumns; |
| 2932 | if (smallIterations < 100) |
| 2933 | smallNumberColumns *= 1.2; |
| 2934 | else |
| 2935 | smallNumberColumns *= 1.1; |
| 2936 | smallNumberColumns = CoinMin(smallNumberColumns, numberColumns); |
| 2937 | if (smallIterations < 200 && smallNumberColumns > 2 * saveSmallNumber) { |
| 2938 | // try kicking it |
| 2939 | emergencyMode = true; |
| 2940 | smallNumberColumns = numberColumns; |
| 2941 | } |
| 2942 | // smallNumberColumns = CoinMin(smallNumberColumns, 3*saveSmallNumber); |
| 2943 | char line[100]; |
| 2944 | sprintf(line, "sample size increased from %d to %d", |
| 2945 | oldNumber, smallNumberColumns); |
| 2946 | handler_->message(CLP_GENERAL, messages_) |
| 2947 | << line |
| 2948 | << CoinMessageEol; |
| 2949 | } |
| 2950 | // move solution back |
| 2951 | const double *solution = small.primalColumnSolution(); |
| 2952 | for (iColumn = 0; iColumn < numberSort; iColumn++) { |
| 2953 | int kColumn = sort[iColumn]; |
| 2954 | model2->setColumnStatus(kColumn, small.getColumnStatus(iColumn)); |
| 2955 | fullSolution[kColumn] = solution[iColumn]; |
| 2956 | } |
| 2957 | for (iRow = 0; iRow < numberRows; iRow++) |
| 2958 | model2->setRowStatus(iRow, small.getRowStatus(iRow)); |
| 2959 | CoinMemcpyN(small.primalRowSolution(), |
| 2960 | numberRows, model2->primalRowSolution()); |
| 2961 | double sumArtificials = 0.0; |
| 2962 | for (i = 0; i < numberArtificials; i++) |
| 2963 | sumArtificials += fullSolution[i + originalNumberColumns]; |
| 2964 | if (sumArtificials && iPass > 5 && sumArtificials >= lastSumArtificials) { |
| 2965 | // increase costs |
| 2966 | double *cost = model2->objective() + originalNumberColumns; |
| 2967 | double newCost = CoinMin(1.0e10, cost[0] * 1.5); |
| 2968 | for (i = 0; i < numberArtificials; i++) |
| 2969 | cost[i] = newCost; |
| 2970 | } |
| 2971 | lastSumArtificials = sumArtificials; |
| 2972 | // get reduced cost for large problem |
| 2973 | double *djs = model2->dualColumnSolution(); |
| 2974 | CoinMemcpyN(model2->objective(), numberColumns, djs); |
| 2975 | model2->clpMatrix()->transposeTimes(-1.0, small.dualRowSolution(), djs); |
| 2976 | int numberNegative = 0; |
| 2977 | double sumNegative = 0.0; |
| 2978 | // now massage weight so all basic in plus good djs |
| 2979 | // first count and do basic |
| 2980 | numberSort = 0; |
| 2981 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
| 2982 | double dj = djs[iColumn] * optimizationDirection_; |
| 2983 | double value = fullSolution[iColumn]; |
| 2984 | if (model2->getColumnStatus(iColumn) == ClpSimplex::basic) { |
| 2985 | sort[numberSort++] = iColumn; |
| 2986 | } else if (dj < -dualTolerance_ && value < columnUpper[iColumn]) { |
| 2987 | numberNegative++; |
| 2988 | sumNegative -= dj; |
| 2989 | } else if (dj > dualTolerance_ && value > columnLower[iColumn]) { |
| 2990 | numberNegative++; |
| 2991 | sumNegative += dj; |
| 2992 | } |
| 2993 | } |
| 2994 | handler_->message(CLP_SPRINT, messages_) |
| 2995 | << iPass + 1 << small.numberIterations() << small.objectiveValue() << sumNegative |
| 2996 | << numberNegative |
| 2997 | << CoinMessageEol; |
| 2998 | if (sumArtificials < 1.0e-8 && originalMaxSprintPass >= 0) { |
| 2999 | maxSprintPass = iPass + originalMaxSprintPass; |
| 3000 | originalMaxSprintPass = -1; |
| 3001 | } |
| 3002 | if (iPass > 20) |
| 3003 | sumArtificials = 0.0; |
| 3004 | if ((small.objectiveValue() * optimizationDirection_ > lastObjective[1] - 1.0e-7 && iPass > 15 && sumArtificials < 1.0e-8 && maxSprintPass < 200) || (!small.numberIterations() && iPass) || iPass == maxSprintPass - 1 || small.status() == 3) { |
3026 | | break; // finished |
3027 | | } else { |
3028 | | lastObjective[1] = lastObjective[0]; |
3029 | | lastObjective[0] = small.objectiveValue() * optimizationDirection_; |
3030 | | double tolerance; |
3031 | | double averageNegDj = sumNegative / static_cast<double> (numberNegative + 1); |
3032 | | if (numberNegative + numberSort > smallNumberColumns && false) |
3033 | | tolerance = -dualTolerance_; |
3034 | | else |
3035 | | tolerance = 10.0 * averageNegDj; |
3036 | | if (emergencyMode) |
3037 | | tolerance=1.0e100; |
3038 | | int saveN = numberSort; |
3039 | | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
3040 | | double dj = djs[iColumn] * optimizationDirection_; |
3041 | | double value = fullSolution[iColumn]; |
3042 | | if (model2->getColumnStatus(iColumn) != ClpSimplex::basic) { |
3043 | | if (dj < -dualTolerance_ && value < columnUpper[iColumn]) |
3044 | | dj = dj; |
3045 | | else if (dj > dualTolerance_ && value > columnLower[iColumn]) |
3046 | | dj = -dj; |
3047 | | else if (columnUpper[iColumn] > columnLower[iColumn]) |
3048 | | dj = fabs(dj); |
3049 | | else |
3050 | | dj = 1.0e50; |
3051 | | if (dj < tolerance) { |
3052 | | weight[numberSort] = dj; |
3053 | | sort[numberSort++] = iColumn; |
3054 | | } |
3055 | | } |
3056 | | } |
3057 | | // sort |
3058 | | CoinSort_2(weight + saveN, weight + numberSort, sort + saveN); |
3059 | | if (numberSort<smallNumberColumns) |
3060 | | printf("using %d columns not %d\n",numberSort,smallNumberColumns); |
3061 | | numberSort = CoinMin(smallNumberColumns, numberSort); |
3062 | | // try singletons |
3063 | | char * markX = new char[numberColumns]; |
3064 | | memset(markX,0,numberColumns); |
3065 | | for (int i=0;i<numberSort;i++) |
3066 | | markX[sort[i]]=1; |
3067 | | int n=numberSort; |
3068 | | for (int i=0;i<numberColumns;i++) { |
3069 | | if (columnLength[i]==1&&!markX[i]) |
3070 | | sort[numberSort++]=i; |
3071 | | } |
3072 | | if (n<numberSort) |
3073 | | printf("%d slacks added\n",numberSort-n); |
3074 | | delete [] markX; |
3075 | | } |
| 3006 | break; // finished |
| 3007 | } else { |
| 3008 | lastObjective[1] = lastObjective[0]; |
| 3009 | lastObjective[0] = small.objectiveValue() * optimizationDirection_; |
| 3010 | double tolerance; |
| 3011 | double averageNegDj = sumNegative / static_cast< double >(numberNegative + 1); |
| 3012 | if (numberNegative + numberSort > smallNumberColumns && false) |
| 3013 | tolerance = -dualTolerance_; |
| 3014 | else |
| 3015 | tolerance = 10.0 * averageNegDj; |
| 3016 | if (emergencyMode) |
| 3017 | tolerance = 1.0e100; |
| 3018 | int saveN = numberSort; |
| 3019 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
| 3020 | double dj = djs[iColumn] * optimizationDirection_; |
| 3021 | double value = fullSolution[iColumn]; |
| 3022 | if (model2->getColumnStatus(iColumn) != ClpSimplex::basic) { |
| 3023 | if (dj < -dualTolerance_ && value < columnUpper[iColumn]) |
| 3024 | dj = dj; |
| 3025 | else if (dj > dualTolerance_ && value > columnLower[iColumn]) |
| 3026 | dj = -dj; |
| 3027 | else if (columnUpper[iColumn] > columnLower[iColumn]) |
| 3028 | dj = fabs(dj); |
| 3029 | else |
| 3030 | dj = 1.0e50; |
| 3031 | if (dj < tolerance) { |
| 3032 | weight[numberSort] = dj; |
| 3033 | sort[numberSort++] = iColumn; |
| 3034 | } |