1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the HiGHS linear optimization suite */
4 /* */
5 /* Written and engineered 2008-2020 at the University of Edinburgh */
6 /* */
7 /* Available as open-source under the MIT License */
8 /* */
9 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
10 /**@file simplex/HighsSimplexAnalysis.cpp
11 * @brief
12 * @author Julian Hall, Ivet Galabova, Qi Huangfu and Michael Feldmeier
13 */
14 #include <cmath>
15 //#include <cstdio>
16 #include "HConfig.h"
17 #include "simplex/FactorTimer.h"
18 #include "simplex/HFactor.h"
19 #include "simplex/HighsSimplexAnalysis.h"
20 #include "simplex/SimplexTimer.h"
21
setup(const HighsLp & lp,const HighsOptions & options,const int simplex_iteration_count_)22 void HighsSimplexAnalysis::setup(const HighsLp& lp, const HighsOptions& options,
23 const int simplex_iteration_count_) {
24 // Copy Problem size
25 numRow = lp.numRow_;
26 numCol = lp.numCol_;
27 numTot = numRow + numCol;
28 // Copy tolerances from options
29 allow_dual_steepest_edge_to_devex_switch =
30 options.simplex_dual_edge_weight_strategy ==
31 SIMPLEX_DUAL_EDGE_WEIGHT_STRATEGY_CHOOSE;
32 dual_steepest_edge_weight_log_error_threshold =
33 options.dual_steepest_edge_weight_log_error_threshold;
34 //
35 AnIterIt0 = simplex_iteration_count_;
36 AnIterCostlyDseFq = 0;
37 AnIterNumCostlyDseIt = 0;
38 // Copy messaging parameter from options
39 messaging(options.logfile, options.output, options.message_level);
40 // Initialise the densities
41 col_aq_density = 0;
42 row_ep_density = 0;
43 row_ap_density = 0;
44 row_DSE_density = 0;
45 col_BFRT_density = 0;
46 primal_col_density = 0;
47 // Set the row_dual_density to 1 since it's assumed all costs are at
48 // least perturbed from zero, if not initially nonzero
49 dual_col_density = 1;
50 // Set up the data structures for scatter data
51 tran_stage.resize(NUM_TRAN_STAGE_TYPE);
52 tran_stage[TRAN_STAGE_FTRAN_LOWER].name_ = "FTRAN lower";
53 tran_stage[TRAN_STAGE_FTRAN_UPPER_FT].name_ = "FTRAN upper FT";
54 tran_stage[TRAN_STAGE_FTRAN_UPPER].name_ = "FTRAN upper";
55 tran_stage[TRAN_STAGE_BTRAN_UPPER].name_ = "BTRAN upper";
56 tran_stage[TRAN_STAGE_BTRAN_UPPER_FT].name_ = "BTRAN upper FT";
57 tran_stage[TRAN_STAGE_BTRAN_LOWER].name_ = "BTRAN lower";
58 for (int tran_stage_type = 0; tran_stage_type < NUM_TRAN_STAGE_TYPE;
59 tran_stage_type++) {
60 TranStageAnalysis& stage = tran_stage[tran_stage_type];
61 initialiseScatterData(20, stage.rhs_density_);
62 stage.num_decision_ = 0;
63 stage.num_wrong_original_sparse_decision_ = 0;
64 stage.num_wrong_original_hyper_decision_ = 0;
65 stage.num_wrong_new_sparse_decision_ = 0;
66 stage.num_wrong_new_hyper_decision_ = 0;
67 }
68 original_start_density_tolerance.resize(NUM_TRAN_STAGE_TYPE);
69 new_start_density_tolerance.resize(NUM_TRAN_STAGE_TYPE);
70 historical_density_tolerance.resize(NUM_TRAN_STAGE_TYPE);
71 predicted_density_tolerance.resize(NUM_TRAN_STAGE_TYPE);
72
73 for (int tran_stage_type = 0; tran_stage_type < NUM_TRAN_STAGE_TYPE;
74 tran_stage_type++) {
75 original_start_density_tolerance[tran_stage_type] = 0.05;
76 new_start_density_tolerance[tran_stage_type] = 0.05;
77 }
78 historical_density_tolerance[TRAN_STAGE_FTRAN_LOWER] = 0.15;
79 historical_density_tolerance[TRAN_STAGE_FTRAN_UPPER] = 0.10;
80 historical_density_tolerance[TRAN_STAGE_BTRAN_UPPER] = 0.10;
81 historical_density_tolerance[TRAN_STAGE_BTRAN_LOWER] = 0.15;
82 predicted_density_tolerance[TRAN_STAGE_FTRAN_LOWER] = 0.10;
83 predicted_density_tolerance[TRAN_STAGE_FTRAN_UPPER] = 0.10;
84 predicted_density_tolerance[TRAN_STAGE_BTRAN_UPPER] = 0.10;
85 predicted_density_tolerance[TRAN_STAGE_BTRAN_LOWER] = 0.10;
86
87 // Initialise the measures used to analyse accuracy of steepest edge weights
88 //
89 const int dual_edge_weight_strategy =
90 options.simplex_dual_edge_weight_strategy;
91 if (dual_edge_weight_strategy == SIMPLEX_DUAL_EDGE_WEIGHT_STRATEGY_CHOOSE ||
92 dual_edge_weight_strategy ==
93 SIMPLEX_DUAL_EDGE_WEIGHT_STRATEGY_STEEPEST_EDGE ||
94 dual_edge_weight_strategy ==
95 SIMPLEX_DUAL_EDGE_WEIGHT_STRATEGY_STEEPEST_EDGE_UNIT_INITIAL) {
96 // Initialise the measures used to analyse accuracy of steepest edge weights
97 num_dual_steepest_edge_weight_check = 0;
98 num_dual_steepest_edge_weight_reject = 0;
99 num_wrong_low_dual_steepest_edge_weight = 0;
100 num_wrong_high_dual_steepest_edge_weight = 0;
101 average_frequency_low_dual_steepest_edge_weight = 0;
102 average_frequency_high_dual_steepest_edge_weight = 0;
103 average_log_low_dual_steepest_edge_weight_error = 0;
104 average_log_high_dual_steepest_edge_weight_error = 0;
105 max_average_frequency_low_dual_steepest_edge_weight = 0;
106 max_average_frequency_high_dual_steepest_edge_weight = 0;
107 max_sum_average_frequency_extreme_dual_steepest_edge_weight = 0;
108 max_average_log_low_dual_steepest_edge_weight_error = 0;
109 max_average_log_high_dual_steepest_edge_weight_error = 0;
110 max_sum_average_log_extreme_dual_steepest_edge_weight_error = 0;
111 }
112 num_devex_framework = 0;
113
114 num_iteration_report_since_last_header = -1;
115 num_invert_report_since_last_header = -1;
116
117 // Set following averages to illegal values so that first average is
118 // set equal to first value
119 average_num_threads = -1;
120 average_fraction_of_possible_minor_iterations_performed = -1;
121 sum_multi_chosen = 0;
122 sum_multi_finished = 0;
123
124 #ifdef HiGHSDEV
125 AnIterPrevIt = simplex_iteration_count_;
126
127 SimplexTimer simplex_timer;
128 for (HighsTimerClock& clock : thread_simplex_clocks)
129 simplex_timer.initialiseSimplexClocks(clock);
130
131 FactorTimer factor_timer;
132 for (HighsTimerClock& clock : thread_factor_clocks)
133 factor_timer.initialiseFactorClocks(clock);
134
135 AnIterOpRec* AnIter;
136 AnIter = &AnIterOp[ANALYSIS_OPERATION_TYPE_BTRAN_EP];
137 AnIter->AnIterOpName = "BTRAN e_p";
138 AnIter = &AnIterOp[ANALYSIS_OPERATION_TYPE_PRICE_AP];
139 AnIter->AnIterOpName = "PRICE a_p";
140 AnIter = &AnIterOp[ANALYSIS_OPERATION_TYPE_BTRAN_FULL];
141 AnIter->AnIterOpName = "BTRAN Full";
142 AnIter = &AnIterOp[ANALYSIS_OPERATION_TYPE_PRICE_FULL];
143 AnIter->AnIterOpName = "PRICE Full";
144 AnIter = &AnIterOp[ANALYSIS_OPERATION_TYPE_FTRAN];
145 AnIter->AnIterOpName = "FTRAN";
146 AnIter = &AnIterOp[ANALYSIS_OPERATION_TYPE_FTRAN_BFRT];
147 AnIter->AnIterOpName = "FTRAN BFRT";
148 AnIter = &AnIterOp[ANALYSIS_OPERATION_TYPE_FTRAN_DSE];
149 AnIter->AnIterOpName = "FTRAN DSE";
150 for (int k = 0; k < NUM_ANALYSIS_OPERATION_TYPE; k++) {
151 AnIter = &AnIterOp[k];
152 if ((k == ANALYSIS_OPERATION_TYPE_PRICE_AP) ||
153 (k == ANALYSIS_OPERATION_TYPE_PRICE_FULL)) {
154 AnIter->AnIterOpHyperCANCEL = 1.0;
155 AnIter->AnIterOpHyperTRAN = 1.0;
156 AnIter->AnIterOpRsDim = numCol;
157 } else {
158 if ((k == ANALYSIS_OPERATION_TYPE_BTRAN_EP) ||
159 (k == ANALYSIS_OPERATION_TYPE_BTRAN_FULL)) {
160 AnIter->AnIterOpHyperCANCEL = hyperCANCEL;
161 AnIter->AnIterOpHyperTRAN = hyperBTRANU;
162 } else {
163 AnIter->AnIterOpHyperCANCEL = hyperCANCEL;
164 AnIter->AnIterOpHyperTRAN = hyperFTRANL;
165 }
166 AnIter->AnIterOpRsDim = numRow;
167 }
168 AnIter->AnIterOpNumCa = 0;
169 AnIter->AnIterOpNumHyperOp = 0;
170 AnIter->AnIterOpNumHyperRs = 0;
171 AnIter->AnIterOpSumLog10RsDensity = 0;
172 initialiseValueDistribution("", "density ", 1e-8, 1.0, 10.0,
173 AnIter->AnIterOp_density);
174 }
175 int last_invert_hint = INVERT_HINT_Count - 1;
176 for (int k = 1; k <= last_invert_hint; k++) AnIterNumInvert[k] = 0;
177 num_col_price = 0;
178 num_row_price = 0;
179 num_row_price_with_switch = 0;
180 int last_dual_edge_weight_mode = (int)DualEdgeWeightMode::STEEPEST_EDGE;
181 for (int k = 0; k <= last_dual_edge_weight_mode; k++) AnIterNumEdWtIt[k] = 0;
182 AnIterTraceNumRec = 0;
183 AnIterTraceIterDl = 1;
184 AnIterTraceRec* lcAnIter = &AnIterTrace[0];
185 lcAnIter->AnIterTraceIter = AnIterIt0;
186 lcAnIter->AnIterTraceTime = timer_->getWallTime();
187 initialiseValueDistribution("Primal step summary", "", 1e-16, 1e16, 10.0,
188 primal_step_distribution);
189 initialiseValueDistribution("Dual step summary", "", 1e-16, 1e16, 10.0,
190 dual_step_distribution);
191 initialiseValueDistribution("Simplex pivot summary", "", 1e-8, 1e16, 10.0,
192 simplex_pivot_distribution);
193 initialiseValueDistribution("Factor pivot threshold summary", "",
194 min_pivot_threshold, max_pivot_threshold,
195 pivot_threshold_change_factor,
196 factor_pivot_threshold_distribution);
197 initialiseValueDistribution("Numerical trouble summary", "", 1e-16, 1.0, 10.0,
198 numerical_trouble_distribution);
199 initialiseValueDistribution("", "1 ", 1e-16, 1e16, 10.0,
200 cost_perturbation1_distribution);
201 initialiseValueDistribution("", "2 ", 1e-16, 1e16, 10.0,
202 cost_perturbation2_distribution);
203 initialiseValueDistribution("FTRAN upper sparse summary - before", "", 1e-8,
204 1.0, 10.0, before_ftran_upper_sparse_density);
205 initialiseValueDistribution("FTRAN upper sparse summary - after", "", 1e-8,
206 1.0, 10.0, before_ftran_upper_hyper_density);
207 initialiseValueDistribution("FTRAN upper hyper-sparse summary - before", "",
208 1e-8, 1.0, 10.0, ftran_upper_sparse_density);
209 initialiseValueDistribution("FTRAN upper hyper-sparse summary - after", "",
210 1e-8, 1.0, 10.0, ftran_upper_hyper_density);
211 initialiseValueDistribution("Cleanup dual change summary", "", 1e-16, 1e16,
212 10.0, cleanup_dual_change_distribution);
213 initialiseValueDistribution("Cleanup primal change summary", "", 1e-16, 1e16,
214 10.0, cleanup_primal_step_distribution);
215 initialiseValueDistribution("Cleanup primal step summary", "", 1e-16, 1e16,
216 10.0, cleanup_primal_change_distribution);
217 initialiseValueDistribution("Cleanup dual step summary", "", 1e-16, 1e16,
218 10.0, cleanup_dual_step_distribution);
219 #endif
220 }
221
messaging(FILE * logfile_,FILE * output_,const int message_level_)222 void HighsSimplexAnalysis::messaging(FILE* logfile_, FILE* output_,
223 const int message_level_) {
224 logfile = logfile_;
225 output = output_;
226 message_level = message_level_;
227 }
228
updateOperationResultDensity(const double local_density,double & density)229 void HighsSimplexAnalysis::updateOperationResultDensity(
230 const double local_density, double& density) {
231 density = (1 - running_average_multiplier) * density +
232 running_average_multiplier * local_density;
233 }
234
iterationReport()235 void HighsSimplexAnalysis::iterationReport() {
236 if (!(iteration_report_message_level & message_level)) return;
237 const bool header = (num_iteration_report_since_last_header < 0) ||
238 (num_iteration_report_since_last_header > 49);
239 if (header) {
240 iterationReport(header);
241 num_iteration_report_since_last_header = 0;
242 }
243 iterationReport(false);
244 }
245
invertReport()246 void HighsSimplexAnalysis::invertReport() {
247 if (!(invert_report_message_level & message_level)) return;
248 const bool header = (num_invert_report_since_last_header < 0) ||
249 (num_invert_report_since_last_header > 49) ||
250 (num_iteration_report_since_last_header >= 0);
251 if (header) {
252 invertReport(header);
253 num_invert_report_since_last_header = 0;
254 }
255 invertReport(false);
256 // Force an iteration report header if this is an INVERT report without an
257 // invert_hint
258 if (!invert_hint) num_iteration_report_since_last_header = -1;
259 }
260
invertReport(const bool header)261 void HighsSimplexAnalysis::invertReport(const bool header) {
262 if (!(invert_report_message_level & message_level)) return;
263 reportAlgorithmPhaseIterationObjective(header, invert_report_message_level);
264 #ifdef HiGHSDEV
265 if (simplex_strategy == SIMPLEX_STRATEGY_DUAL_MULTI) {
266 // Report on threads and PAMI
267 reportThreads(header, invert_report_message_level);
268 reportMulti(header, invert_report_message_level);
269 }
270 reportDensity(header, invert_report_message_level);
271 reportInvert(header, invert_report_message_level);
272 // reportCondition(header, invert_report_message_level);
273 #endif
274 reportInfeasibility(header, invert_report_message_level);
275 HighsPrintMessage(output, message_level, invert_report_message_level, "\n");
276 if (!header) num_invert_report_since_last_header++;
277 }
278
dualSteepestEdgeWeightError(const double computed_edge_weight,const double updated_edge_weight)279 void HighsSimplexAnalysis::dualSteepestEdgeWeightError(
280 const double computed_edge_weight, const double updated_edge_weight) {
281 const bool accept_weight =
282 updated_edge_weight >= accept_weight_threshold * computed_edge_weight;
283 int low_weight_error = 0;
284 int high_weight_error = 0;
285 double weight_error;
286 #ifdef HiGHSDEV
287 string error_type = " OK";
288 #endif
289 num_dual_steepest_edge_weight_check++;
290 if (!accept_weight) num_dual_steepest_edge_weight_reject++;
291 if (updated_edge_weight < computed_edge_weight) {
292 // Updated weight is low
293 weight_error = computed_edge_weight / updated_edge_weight;
294 if (weight_error > weight_error_threshold) {
295 low_weight_error = 1;
296 #ifdef HiGHSDEV
297 error_type = " Low";
298 #endif
299 }
300 average_log_low_dual_steepest_edge_weight_error =
301 0.99 * average_log_low_dual_steepest_edge_weight_error +
302 0.01 * log(weight_error);
303 } else {
304 // Updated weight is correct or high
305 weight_error = updated_edge_weight / computed_edge_weight;
306 if (weight_error > weight_error_threshold) {
307 high_weight_error = 1;
308 #ifdef HiGHSDEV
309 error_type = "High";
310 #endif
311 }
312 average_log_high_dual_steepest_edge_weight_error =
313 0.99 * average_log_high_dual_steepest_edge_weight_error +
314 0.01 * log(weight_error);
315 }
316 average_frequency_low_dual_steepest_edge_weight =
317 0.99 * average_frequency_low_dual_steepest_edge_weight +
318 0.01 * low_weight_error;
319 average_frequency_high_dual_steepest_edge_weight =
320 0.99 * average_frequency_high_dual_steepest_edge_weight +
321 0.01 * high_weight_error;
322 max_average_frequency_low_dual_steepest_edge_weight =
323 max(max_average_frequency_low_dual_steepest_edge_weight,
324 average_frequency_low_dual_steepest_edge_weight);
325 max_average_frequency_high_dual_steepest_edge_weight =
326 max(max_average_frequency_high_dual_steepest_edge_weight,
327 average_frequency_high_dual_steepest_edge_weight);
328 max_sum_average_frequency_extreme_dual_steepest_edge_weight =
329 max(max_sum_average_frequency_extreme_dual_steepest_edge_weight,
330 average_frequency_low_dual_steepest_edge_weight +
331 average_frequency_high_dual_steepest_edge_weight);
332 max_average_log_low_dual_steepest_edge_weight_error =
333 max(max_average_log_low_dual_steepest_edge_weight_error,
334 average_log_low_dual_steepest_edge_weight_error);
335 max_average_log_high_dual_steepest_edge_weight_error =
336 max(max_average_log_high_dual_steepest_edge_weight_error,
337 average_log_high_dual_steepest_edge_weight_error);
338 max_sum_average_log_extreme_dual_steepest_edge_weight_error =
339 max(max_sum_average_log_extreme_dual_steepest_edge_weight_error,
340 average_log_low_dual_steepest_edge_weight_error +
341 average_log_high_dual_steepest_edge_weight_error);
342 #ifdef HiGHSDEV
343 const bool report_weight_error = false;
344 if (report_weight_error && weight_error > 0.5 * weight_error_threshold) {
345 printf(
346 "DSE Wt Ck |%8d| OK = %1d (%4d / %6d) (c %10.4g, u %10.4g, er %10.4g - "
347 "%s): Low (Fq %10.4g, Er %10.4g); High (Fq%10.4g, Er%10.4g) | %10.4g "
348 "%10.4g %10.4g %10.4g %10.4g %10.4g\n",
349 simplex_iteration_count, accept_weight,
350 num_dual_steepest_edge_weight_check,
351 num_dual_steepest_edge_weight_reject, computed_edge_weight,
352 updated_edge_weight, weight_error, error_type.c_str(),
353 average_frequency_low_dual_steepest_edge_weight,
354 average_log_low_dual_steepest_edge_weight_error,
355 average_frequency_high_dual_steepest_edge_weight,
356 average_log_high_dual_steepest_edge_weight_error,
357 max_average_frequency_low_dual_steepest_edge_weight,
358 max_average_frequency_high_dual_steepest_edge_weight,
359 max_sum_average_frequency_extreme_dual_steepest_edge_weight,
360 max_average_log_low_dual_steepest_edge_weight_error,
361 max_average_log_high_dual_steepest_edge_weight_error,
362 max_sum_average_log_extreme_dual_steepest_edge_weight_error);
363 }
364 #endif
365 }
366
switchToDevex()367 bool HighsSimplexAnalysis::switchToDevex() {
368 bool switch_to_devex = false;
369 // Firstly consider switching on the basis of NLA cost
370 double AnIterCostlyDseMeasureDen;
371 AnIterCostlyDseMeasureDen =
372 max(max(row_ep_density, col_aq_density), row_ap_density);
373 if (AnIterCostlyDseMeasureDen > 0) {
374 AnIterCostlyDseMeasure = row_DSE_density / AnIterCostlyDseMeasureDen;
375 AnIterCostlyDseMeasure = AnIterCostlyDseMeasure * AnIterCostlyDseMeasure;
376 } else {
377 AnIterCostlyDseMeasure = 0;
378 }
379 bool CostlyDseIt = AnIterCostlyDseMeasure > AnIterCostlyDseMeasureLimit &&
380 row_DSE_density > AnIterCostlyDseMnDensity;
381 AnIterCostlyDseFq = (1 - running_average_multiplier) * AnIterCostlyDseFq;
382 if (CostlyDseIt) {
383 AnIterNumCostlyDseIt++;
384 AnIterCostlyDseFq += running_average_multiplier * 1.0;
385 int lcNumIter = simplex_iteration_count - AnIterIt0;
386 // Switch to Devex if at least 5% of the (at least) 0.1NumTot iterations
387 // have been costly
388 switch_to_devex =
389 allow_dual_steepest_edge_to_devex_switch &&
390 (AnIterNumCostlyDseIt > lcNumIter * AnIterFracNumCostlyDseItbfSw) &&
391 (lcNumIter > AnIterFracNumTot_ItBfSw * numTot);
392 #ifdef HiGHSDEV
393 if (switch_to_devex) {
394 HighsLogMessage(
395 logfile, HighsMessageType::INFO,
396 "Switch from DSE to Devex after %d costly DSE iterations of %d: "
397 "C_Aq_Density = %11.4g; R_Ep_Density = %11.4g; R_Ap_Density = "
398 "%11.4g; DSE_Density = %11.4g",
399 AnIterNumCostlyDseIt, lcNumIter, col_aq_density, row_ep_density,
400 row_ap_density, row_DSE_density);
401 }
402 #endif
403 }
404 if (!switch_to_devex) {
405 // Secondly consider switching on the basis of weight accuracy
406 double dse_weight_error_measure =
407 average_log_low_dual_steepest_edge_weight_error +
408 average_log_high_dual_steepest_edge_weight_error;
409 double dse_weight_error_threshold =
410 dual_steepest_edge_weight_log_error_threshold;
411 switch_to_devex = allow_dual_steepest_edge_to_devex_switch &&
412 dse_weight_error_measure > dse_weight_error_threshold;
413 #ifdef HiGHSDEV
414 if (switch_to_devex) {
415 HighsLogMessage(logfile, HighsMessageType::INFO,
416 "Switch from DSE to Devex with log error measure of %g > "
417 "%g = threshold",
418 dse_weight_error_measure, dse_weight_error_threshold);
419 }
420 #endif
421 }
422 return switch_to_devex;
423 }
424
predictEndDensity(const int tran_stage_type,const double start_density,double & end_density)425 bool HighsSimplexAnalysis::predictEndDensity(const int tran_stage_type,
426 const double start_density,
427 double& end_density) {
428 return predictFromScatterData(tran_stage[tran_stage_type].rhs_density_,
429 start_density, end_density);
430 }
431
afterTranStage(const int tran_stage_type,const double start_density,const double end_density,const double historical_density,const double predicted_end_density,const bool use_solve_sparse_original_HFactor_logic,const bool use_solve_sparse_new_HFactor_logic)432 void HighsSimplexAnalysis::afterTranStage(
433 const int tran_stage_type, const double start_density,
434 const double end_density, const double historical_density,
435 const double predicted_end_density,
436 const bool use_solve_sparse_original_HFactor_logic,
437 const bool use_solve_sparse_new_HFactor_logic) {
438 TranStageAnalysis& stage = tran_stage[tran_stage_type];
439 const double rp = false;
440 if (predicted_end_density > 0) {
441 stage.num_decision_++;
442 if (end_density <= max_hyper_density) {
443 // Should have done hyper-sparse TRAN
444 if (use_solve_sparse_original_HFactor_logic) {
445 // Original logic makes wrong decision to use sparse TRAN
446 if (rp) {
447 printf("Original: Wrong sparse: ");
448 const double start_density_tolerance =
449 original_start_density_tolerance[tran_stage_type];
450 const double this_historical_density_tolerance =
451 historical_density_tolerance[tran_stage_type];
452 if (start_density > start_density_tolerance) {
453 printf("(start = %10.4g > %4.2f) or ", start_density,
454 start_density_tolerance);
455 } else {
456 printf(" start = %10.4g ", start_density);
457 }
458 if (historical_density > this_historical_density_tolerance) {
459 printf("(historical = %10.4g > %4.2f); ", historical_density,
460 this_historical_density_tolerance);
461 } else {
462 printf(" historical = %10.4g ", historical_density);
463 }
464 printf("end = %10.4g", end_density);
465 if (end_density < 0.1 * historical_density) printf(" !! OG");
466 printf("\n");
467 }
468 stage.num_wrong_original_sparse_decision_++;
469 }
470 if (use_solve_sparse_new_HFactor_logic) {
471 // New logic makes wrong decision to use sparse TRAN
472 if (rp) {
473 printf("New : Wrong sparse: ");
474 const double start_density_tolerance =
475 original_start_density_tolerance[tran_stage_type];
476 const double end_density_tolerance =
477 predicted_density_tolerance[tran_stage_type];
478 if (start_density > start_density_tolerance) {
479 printf("(start = %10.4g > %4.2f) or ", start_density,
480 start_density_tolerance);
481 } else {
482 printf(" start = %10.4g ", start_density);
483 }
484 if (predicted_end_density > end_density_tolerance) {
485 printf("( predicted = %10.4g > %4.2f); ", predicted_end_density,
486 end_density_tolerance);
487 } else {
488 printf(" predicted = %10.4g ", predicted_end_density);
489 }
490 printf("end = %10.4g", end_density);
491 if (end_density < 0.1 * predicted_end_density) printf(" !! NW");
492 printf("\n");
493 }
494 stage.num_wrong_new_sparse_decision_++;
495 }
496 } else {
497 // Should have done sparse TRAN
498 if (!use_solve_sparse_original_HFactor_logic) {
499 // Original logic makes wrong decision to use hyper TRAN
500 if (rp) {
501 printf(
502 "Original: Wrong hyper: (start = %10.4g <= %4.2f) and "
503 "(historical = %10.4g <= %4.2f); end = %10.4g",
504 start_density, original_start_density_tolerance[tran_stage_type],
505 historical_density, historical_density_tolerance[tran_stage_type],
506 end_density);
507 if (end_density > 10.0 * historical_density) printf(" !! OG");
508 printf("\n");
509 }
510 stage.num_wrong_original_hyper_decision_++;
511 }
512 if (!use_solve_sparse_new_HFactor_logic) {
513 // New logic makes wrong decision to use hyper TRAN
514 if (rp) {
515 printf(
516 "New : Wrong hyper: (start = %10.4g <= %4.2f) and ( "
517 "predicted = %10.4g <= %4.2f); end = %10.4g",
518 start_density, new_start_density_tolerance[tran_stage_type],
519 predicted_end_density,
520 predicted_density_tolerance[tran_stage_type], end_density);
521 if (end_density > 10.0 * predicted_end_density) printf(" !! NW");
522 printf("\n");
523 }
524 stage.num_wrong_new_hyper_decision_++;
525 }
526 }
527 }
528 updateScatterData(start_density, end_density, stage.rhs_density_);
529 regressScatterData(stage.rhs_density_);
530 }
531
summaryReportFactor()532 void HighsSimplexAnalysis::summaryReportFactor() {
533 for (int tran_stage_type = 0; tran_stage_type < NUM_TRAN_STAGE_TYPE;
534 tran_stage_type++) {
535 TranStageAnalysis& stage = tran_stage[tran_stage_type];
536 // printScatterData(stage.name_, stage.rhs_density_);
537 printScatterDataRegressionComparison(stage.name_, stage.rhs_density_);
538 if (!stage.num_decision_) return;
539 printf("Of %10d Sps/Hyper decisions made using regression:\n",
540 stage.num_decision_);
541 printf(
542 " %10d wrong sparseTRAN; %10d wrong hyperTRAN: using original "
543 "logic\n",
544 stage.num_wrong_original_sparse_decision_,
545 stage.num_wrong_original_hyper_decision_);
546 printf(
547 " %10d wrong sparseTRAN; %10d wrong hyperTRAN: using new "
548 "logic\n",
549 stage.num_wrong_new_sparse_decision_,
550 stage.num_wrong_new_hyper_decision_);
551 }
552 }
553
simplexTimerStart(const int simplex_clock,const int thread_id)554 void HighsSimplexAnalysis::simplexTimerStart(const int simplex_clock,
555 const int thread_id) {
556 #ifdef HiGHSDEV
557 thread_simplex_clocks[thread_id].timer_.start(
558 thread_simplex_clocks[thread_id].clock_[simplex_clock]);
559 #endif
560 }
561
simplexTimerStop(const int simplex_clock,const int thread_id)562 void HighsSimplexAnalysis::simplexTimerStop(const int simplex_clock,
563 const int thread_id) {
564 #ifdef HiGHSDEV
565 thread_simplex_clocks[thread_id].timer_.stop(
566 thread_simplex_clocks[thread_id].clock_[simplex_clock]);
567 #endif
568 }
569
simplexTimerRunning(const int simplex_clock,const int thread_id)570 bool HighsSimplexAnalysis::simplexTimerRunning(const int simplex_clock,
571 const int thread_id) {
572 bool argument = false;
573 #ifdef HiGHSDEV
574 argument =
575 thread_simplex_clocks[thread_id]
576 .timer_
577 .clock_start[thread_simplex_clocks[thread_id].clock_[simplex_clock]] <
578 0;
579 #endif
580 return argument;
581 }
582
simplexTimerNumCall(const int simplex_clock,const int thread_id)583 int HighsSimplexAnalysis::simplexTimerNumCall(const int simplex_clock,
584 const int thread_id) {
585 int argument = 0;
586 #ifdef HiGHSDEV
587 argument = thread_simplex_clocks[thread_id].timer_.clock_num_call
588 [thread_simplex_clocks[thread_id].clock_[simplex_clock]];
589 #endif
590 return argument;
591 }
592
simplexTimerRead(const int simplex_clock,const int thread_id)593 double HighsSimplexAnalysis::simplexTimerRead(const int simplex_clock,
594 const int thread_id) {
595 double argument = 0;
596 #ifdef HiGHSDEV
597 argument = thread_simplex_clocks[thread_id].timer_.read(
598 thread_simplex_clocks[thread_id].clock_[simplex_clock]);
599 #endif
600 return argument;
601 }
602
getThreadFactorTimerClockPointer()603 HighsTimerClock* HighsSimplexAnalysis::getThreadFactorTimerClockPointer() {
604 HighsTimerClock* factor_timer_clock_pointer = NULL;
605 #ifdef HiGHSDEV
606 int thread_id = 0;
607 #ifdef OPENMP
608 thread_id = omp_get_thread_num();
609 #endif
610 factor_timer_clock_pointer = &thread_factor_clocks[thread_id];
611 #endif
612 return factor_timer_clock_pointer;
613 }
614
615 #ifdef HiGHSDEV
reportFactorTimer()616 void HighsSimplexAnalysis::reportFactorTimer() {
617 FactorTimer factor_timer;
618 int omp_max_threads = 1;
619 #ifdef OPENMP
620 omp_max_threads = omp_get_max_threads();
621 #endif
622 for (int i = 0; i < omp_max_threads; i++) {
623 // for (HighsTimerClock clock : thread_factor_clocks) {
624 printf("reportFactorTimer: HFactor clocks for OMP thread %d / %d\n", i,
625 omp_max_threads - 1);
626 factor_timer.reportFactorClock(thread_factor_clocks[i]);
627 }
628 if (omp_max_threads > 1) {
629 HighsTimer& timer = thread_factor_clocks[0].timer_;
630 HighsTimerClock all_factor_clocks(timer);
631 vector<int>& clock = all_factor_clocks.clock_;
632 factor_timer.initialiseFactorClocks(all_factor_clocks);
633 for (int i = 0; i < omp_max_threads; i++) {
634 vector<int>& thread_clock = thread_factor_clocks[i].clock_;
635 for (int clock_id = 0; clock_id < FactorNumClock; clock_id++) {
636 int all_factor_iClock = clock[clock_id];
637 int thread_factor_iClock = thread_clock[clock_id];
638 timer.clock_num_call[all_factor_iClock] +=
639 timer.clock_num_call[thread_factor_iClock];
640 timer.clock_time[all_factor_iClock] +=
641 timer.clock_time[thread_factor_iClock];
642 }
643 }
644 printf("reportFactorTimer: HFactor clocks for all %d threads\n",
645 omp_max_threads);
646 factor_timer.reportFactorClock(all_factor_clocks);
647 }
648 }
649
iterationRecord()650 void HighsSimplexAnalysis::iterationRecord() {
651 int AnIterCuIt = simplex_iteration_count;
652 if (invert_hint > 0) AnIterNumInvert[invert_hint]++;
653 if (AnIterCuIt > AnIterPrevIt)
654 AnIterNumEdWtIt[(int)edge_weight_mode] += (AnIterCuIt - AnIterPrevIt);
655
656 AnIterTraceRec& lcAnIter = AnIterTrace[AnIterTraceNumRec];
657 // if (simplex_iteration_count ==
658 // AnIterTraceIterRec[AnIterTraceNumRec]+AnIterTraceIterDl) {
659 if (simplex_iteration_count == lcAnIter.AnIterTraceIter + AnIterTraceIterDl) {
660 if (AnIterTraceNumRec == AN_ITER_TRACE_MX_NUM_REC) {
661 for (int rec = 1; rec <= AN_ITER_TRACE_MX_NUM_REC / 2; rec++)
662 AnIterTrace[rec] = AnIterTrace[2 * rec];
663 AnIterTraceNumRec = AnIterTraceNumRec / 2;
664 AnIterTraceIterDl = AnIterTraceIterDl * 2;
665 } else {
666 AnIterTraceNumRec++;
667 AnIterTraceRec& lcAnIter = AnIterTrace[AnIterTraceNumRec];
668 lcAnIter.AnIterTraceIter = simplex_iteration_count;
669 lcAnIter.AnIterTraceTime = timer_->getWallTime();
670 if (average_fraction_of_possible_minor_iterations_performed > 0) {
671 lcAnIter.AnIterTraceMulti =
672 average_fraction_of_possible_minor_iterations_performed;
673 } else {
674 lcAnIter.AnIterTraceMulti = 0;
675 }
676 lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN] =
677 col_aq_density;
678 lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_BTRAN_EP] =
679 row_ep_density;
680 lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_PRICE_AP] =
681 row_ap_density;
682 lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN_BFRT] =
683 col_aq_density;
684 if (edge_weight_mode == DualEdgeWeightMode::STEEPEST_EDGE) {
685 lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN_DSE] =
686 row_DSE_density;
687 lcAnIter.AnIterTraceCostlyDse = AnIterCostlyDseMeasure;
688 } else {
689 lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN_DSE] = 0;
690 lcAnIter.AnIterTraceCostlyDse = 0;
691 }
692 lcAnIter.AnIterTrace_dual_edge_weight_mode = (int)edge_weight_mode;
693 }
694 }
695 AnIterPrevIt = AnIterCuIt;
696 updateValueDistribution(primal_step, cleanup_primal_step_distribution);
697 updateValueDistribution(dual_step, cleanup_dual_step_distribution);
698 updateValueDistribution(primal_step, primal_step_distribution);
699 updateValueDistribution(dual_step, dual_step_distribution);
700 updateValueDistribution(pivot_value_from_column, simplex_pivot_distribution);
701 updateValueDistribution(factor_pivot_threshold,
702 factor_pivot_threshold_distribution);
703 // Only update the distribution of legal values for
704 // numerical_trouble. Illegal values are set in PAMI since it's not
705 // known in minor iterations
706 if (numerical_trouble >= 0)
707 updateValueDistribution(numerical_trouble, numerical_trouble_distribution);
708 }
709
iterationRecordMajor()710 void HighsSimplexAnalysis::iterationRecordMajor() {
711 sum_multi_chosen += multi_chosen;
712 sum_multi_finished += multi_finished;
713 assert(multi_chosen > 0);
714 const double fraction_of_possible_minor_iterations_performed =
715 1.0 * multi_finished / multi_chosen;
716 if (average_fraction_of_possible_minor_iterations_performed < 0) {
717 average_fraction_of_possible_minor_iterations_performed =
718 fraction_of_possible_minor_iterations_performed;
719 } else {
720 average_fraction_of_possible_minor_iterations_performed =
721 running_average_multiplier *
722 fraction_of_possible_minor_iterations_performed +
723 (1 - running_average_multiplier) *
724 average_fraction_of_possible_minor_iterations_performed;
725 }
726 if (average_num_threads < 0) {
727 average_num_threads = num_threads;
728 } else {
729 average_num_threads =
730 running_average_multiplier * num_threads +
731 (1 - running_average_multiplier) * average_num_threads;
732 }
733 }
734
operationRecordBefore(const int operation_type,const HVector & vector,const double historical_density)735 void HighsSimplexAnalysis::operationRecordBefore(
736 const int operation_type, const HVector& vector,
737 const double historical_density) {
738 operationRecordBefore(operation_type, vector.count, historical_density);
739 }
740
operationRecordBefore(const int operation_type,const int current_count,const double historical_density)741 void HighsSimplexAnalysis::operationRecordBefore(
742 const int operation_type, const int current_count,
743 const double historical_density) {
744 double current_density = 1.0 * current_count / numRow;
745 AnIterOpRec& AnIter = AnIterOp[operation_type];
746 AnIter.AnIterOpNumCa++;
747 if (current_density <= AnIter.AnIterOpHyperCANCEL &&
748 historical_density <= AnIter.AnIterOpHyperTRAN)
749 AnIter.AnIterOpNumHyperOp++;
750 }
751
operationRecordAfter(const int operation_type,const HVector & vector)752 void HighsSimplexAnalysis::operationRecordAfter(const int operation_type,
753 const HVector& vector) {
754 operationRecordAfter(operation_type, vector.count);
755 }
756
operationRecordAfter(const int operation_type,const int result_count)757 void HighsSimplexAnalysis::operationRecordAfter(const int operation_type,
758 const int result_count) {
759 AnIterOpRec& AnIter = AnIterOp[operation_type];
760 const double result_density = 1.0 * result_count / AnIter.AnIterOpRsDim;
761 if (result_density <= hyperRESULT) AnIter.AnIterOpNumHyperRs++;
762 if (result_density > 0) {
763 AnIter.AnIterOpSumLog10RsDensity += log(result_density) / log(10.0);
764 } else {
765 /*
766 // TODO Investigate these zero norms
767 double vectorNorm = 0;
768
769 for (int index = 0; index < AnIter.AnIterOpRsDim; index++) {
770 double vectorValue = vector.array[index];
771 vectorNorm += vectorValue * vectorValue;
772 }
773 vectorNorm = sqrt(vectorNorm);
774 printf("Strange: operation %s has result density = %g: ||vector|| = %g\n",
775 AnIter.AnIterOpName.c_str(), result_density, vectorNorm);
776 */
777 }
778 updateValueDistribution(result_density, AnIter.AnIterOp_density);
779 }
780
summaryReport()781 void HighsSimplexAnalysis::summaryReport() {
782 int AnIterNumIter = simplex_iteration_count - AnIterIt0;
783 if (AnIterNumIter <= 0) return;
784 printf("\nAnalysis of %d iterations (%d to %d)\n", AnIterNumIter,
785 AnIterIt0 + 1, simplex_iteration_count);
786 if (AnIterNumIter <= 0) return;
787 int lc_EdWtNumIter;
788 lc_EdWtNumIter = AnIterNumEdWtIt[(int)DualEdgeWeightMode::STEEPEST_EDGE];
789 if (lc_EdWtNumIter > 0)
790 printf("DSE for %12d (%3d%%) iterations\n", lc_EdWtNumIter,
791 (100 * lc_EdWtNumIter) / AnIterNumIter);
792 lc_EdWtNumIter = AnIterNumEdWtIt[(int)DualEdgeWeightMode::DEVEX];
793 if (lc_EdWtNumIter > 0)
794 printf("Dvx for %12d (%3d%%) iterations\n", lc_EdWtNumIter,
795 (100 * lc_EdWtNumIter) / AnIterNumIter);
796 lc_EdWtNumIter = AnIterNumEdWtIt[(int)DualEdgeWeightMode::DANTZIG];
797 if (lc_EdWtNumIter > 0)
798 printf("Dan for %12d (%3d%%) iterations\n", lc_EdWtNumIter,
799 (100 * lc_EdWtNumIter) / AnIterNumIter);
800 for (int k = 0; k < NUM_ANALYSIS_OPERATION_TYPE; k++) {
801 AnIterOpRec& AnIter = AnIterOp[k];
802 int lcNumCa = AnIter.AnIterOpNumCa;
803 printf("\n%-10s performed %d times\n", AnIter.AnIterOpName.c_str(),
804 AnIter.AnIterOpNumCa);
805 if (lcNumCa > 0) {
806 int lcHyperOp = AnIter.AnIterOpNumHyperOp;
807 int lcHyperRs = AnIter.AnIterOpNumHyperRs;
808 int pctHyperOp = (100 * lcHyperOp) / lcNumCa;
809 int pctHyperRs = (100 * lcHyperRs) / lcNumCa;
810 double lcRsDensity =
811 pow(10.0, AnIter.AnIterOpSumLog10RsDensity / lcNumCa);
812 int lcAnIterOpRsDim = AnIter.AnIterOpRsDim;
813 int lcNumNNz = lcRsDensity * lcAnIterOpRsDim;
814 printf("%12d hyper-sparse operations (%3d%%)\n", lcHyperOp, pctHyperOp);
815 printf("%12d hyper-sparse results (%3d%%)\n", lcHyperRs, pctHyperRs);
816 printf("%12g density of result (%d / %d nonzeros)\n", lcRsDensity,
817 lcNumNNz, lcAnIterOpRsDim);
818 printValueDistribution(AnIter.AnIterOp_density, AnIter.AnIterOpRsDim);
819 }
820 }
821 int NumInvert = 0;
822
823 int last_invert_hint = INVERT_HINT_Count - 1;
824 for (int k = 1; k <= last_invert_hint; k++) NumInvert += AnIterNumInvert[k];
825 if (NumInvert > 0) {
826 int lcNumInvert = 0;
827 printf("\nInvert performed %d times: average frequency = %d\n",
828 NumInvert, AnIterNumIter / NumInvert);
829 lcNumInvert = AnIterNumInvert[INVERT_HINT_UPDATE_LIMIT_REACHED];
830 if (lcNumInvert > 0)
831 printf("%12d (%3d%%) Invert operations due to update limit reached\n",
832 lcNumInvert, (100 * lcNumInvert) / NumInvert);
833 lcNumInvert = AnIterNumInvert[INVERT_HINT_SYNTHETIC_CLOCK_SAYS_INVERT];
834 if (lcNumInvert > 0)
835 printf("%12d (%3d%%) Invert operations due to pseudo-clock\n",
836 lcNumInvert, (100 * lcNumInvert) / NumInvert);
837 lcNumInvert = AnIterNumInvert[INVERT_HINT_POSSIBLY_OPTIMAL];
838 if (lcNumInvert > 0)
839 printf("%12d (%3d%%) Invert operations due to possibly optimal\n",
840 lcNumInvert, (100 * lcNumInvert) / NumInvert);
841 lcNumInvert = AnIterNumInvert[INVERT_HINT_POSSIBLY_PRIMAL_UNBOUNDED];
842 if (lcNumInvert > 0)
843 printf(
844 "%12d (%3d%%) Invert operations due to possibly primal unbounded\n",
845 lcNumInvert, (100 * lcNumInvert) / NumInvert);
846 lcNumInvert = AnIterNumInvert[INVERT_HINT_POSSIBLY_DUAL_UNBOUNDED];
847 if (lcNumInvert > 0)
848 printf("%12d (%3d%%) Invert operations due to possibly dual unbounded\n",
849 lcNumInvert, (100 * lcNumInvert) / NumInvert);
850 lcNumInvert = AnIterNumInvert[INVERT_HINT_POSSIBLY_SINGULAR_BASIS];
851 if (lcNumInvert > 0)
852 printf("%12d (%3d%%) Invert operations due to possibly singular basis\n",
853 lcNumInvert, (100 * lcNumInvert) / NumInvert);
854 lcNumInvert =
855 AnIterNumInvert[INVERT_HINT_PRIMAL_INFEASIBLE_IN_PRIMAL_SIMPLEX];
856 if (lcNumInvert > 0)
857 printf(
858 "%12d (%3d%%) Invert operations due to primal infeasible in primal "
859 "simplex\n",
860 lcNumInvert, (100 * lcNumInvert) / NumInvert);
861 }
862 int suPrice = num_col_price + num_row_price + num_row_price_with_switch;
863 if (suPrice > 0) {
864 printf("\n%12d Price operations:\n", suPrice);
865 printf("%12d Col Price (%3d%%)\n", num_col_price,
866 (100 * num_col_price) / suPrice);
867 printf("%12d Row Price (%3d%%)\n", num_row_price,
868 (100 * num_row_price) / suPrice);
869 printf("%12d Row PriceWSw (%3d%%)\n", num_row_price_with_switch,
870 (100 * num_row_price_with_switch / suPrice));
871 }
872 printf("\n%12d (%3d%%) costly DSE iterations\n", AnIterNumCostlyDseIt,
873 (100 * AnIterNumCostlyDseIt) / AnIterNumIter);
874
875 // Look for any Devex data to summarise
876 if (num_devex_framework) {
877 printf("\nDevex summary\n");
878 printf("%12d Devex frameworks\n", num_devex_framework);
879 printf(
880 "%12d average number of iterations\n",
881 AnIterNumEdWtIt[(int)DualEdgeWeightMode::DEVEX] / num_devex_framework);
882 }
883 // Look for any PAMI data to summarise
884 if (sum_multi_chosen > 0) {
885 const int pct_minor_iterations_performed =
886 (100 * sum_multi_finished) / sum_multi_chosen;
887 printf("\nPAMI summary: for average of %0.1g threads \n",
888 average_num_threads);
889 printf("%12d Major iterations\n", multi_iteration_count);
890 printf("%12d Minor iterations\n", sum_multi_finished);
891 printf(
892 "%12d Total rows chosen: performed %3d%% of possible minor "
893 "iterations\n\n",
894 sum_multi_chosen, pct_minor_iterations_performed);
895 }
896
897 printf("\nCost perturbation summary\n");
898 printValueDistribution(cost_perturbation1_distribution);
899 printValueDistribution(cost_perturbation2_distribution);
900
901 printValueDistribution(before_ftran_upper_sparse_density, numRow);
902 printValueDistribution(ftran_upper_sparse_density, numRow);
903 printValueDistribution(before_ftran_upper_hyper_density, numRow);
904 printValueDistribution(ftran_upper_hyper_density, numRow);
905 printValueDistribution(primal_step_distribution);
906 printValueDistribution(dual_step_distribution);
907 printValueDistribution(simplex_pivot_distribution);
908 printValueDistribution(factor_pivot_threshold_distribution);
909 printValueDistribution(numerical_trouble_distribution);
910 printValueDistribution(cleanup_dual_change_distribution);
911 printValueDistribution(cleanup_primal_step_distribution);
912 printValueDistribution(cleanup_dual_step_distribution);
913 printValueDistribution(cleanup_primal_change_distribution);
914
915 if (AnIterTraceIterDl >= 100) {
916 // Possibly (usually) add a temporary record for the final
917 // iterations: may end up with one more than
918 // AN_ITER_TRACE_MX_NUM_REC records, so ensure that there is
919 // enough space in the arrays
920 //
921 const bool add_extra_record =
922 simplex_iteration_count >
923 AnIterTrace[AnIterTraceNumRec].AnIterTraceIter;
924 if (add_extra_record) {
925 AnIterTraceNumRec++;
926 AnIterTraceRec& lcAnIter = AnIterTrace[AnIterTraceNumRec];
927 lcAnIter.AnIterTraceIter = simplex_iteration_count;
928 lcAnIter.AnIterTraceTime = timer_->getWallTime();
929 if (average_fraction_of_possible_minor_iterations_performed > 0) {
930 lcAnIter.AnIterTraceMulti =
931 average_fraction_of_possible_minor_iterations_performed;
932 } else {
933 lcAnIter.AnIterTraceMulti = 0;
934 }
935 lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN] =
936 col_aq_density;
937 lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_BTRAN_EP] =
938 row_ep_density;
939 lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_PRICE_AP] =
940 row_ap_density;
941 lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN_BFRT] =
942 col_aq_density;
943 if (edge_weight_mode == DualEdgeWeightMode::STEEPEST_EDGE) {
944 lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN_DSE] =
945 row_DSE_density;
946 lcAnIter.AnIterTraceCostlyDse = AnIterCostlyDseMeasure;
947 } else {
948 lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN_DSE] = 0;
949 lcAnIter.AnIterTraceCostlyDse = 0;
950 }
951 lcAnIter.AnIterTrace_dual_edge_weight_mode = (int)edge_weight_mode;
952 }
953 // Determine whether the Multi and steepest edge columns should be reported
954 double su_multi_values = 0;
955 double su_dse_values = 0;
956 for (int rec = 1; rec <= AnIterTraceNumRec; rec++) {
957 AnIterTraceRec& lcAnIter = AnIterTrace[rec];
958 su_multi_values += fabs(lcAnIter.AnIterTraceMulti);
959 su_dse_values +=
960 fabs(lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN_DSE]);
961 }
962 const bool report_multi = su_multi_values > 0;
963 const bool rp_dual_steepest_edge = su_dse_values > 0;
964 printf("\n Iteration speed analysis\n");
965 AnIterTraceRec& lcAnIter = AnIterTrace[0];
966 int fmIter = lcAnIter.AnIterTraceIter;
967 double fmTime = lcAnIter.AnIterTraceTime;
968 printf(" Iter ( FmIter: ToIter) Time Iter/sec ");
969 if (report_multi) printf("| PAMI ");
970 printf("| C_Aq R_Ep R_Ap ");
971 if (rp_dual_steepest_edge) printf(" DSE ");
972 printf("| EdWt ");
973 if (rp_dual_steepest_edge) {
974 printf("| CostlyDse\n");
975 } else {
976 printf("\n");
977 }
978
979 for (int rec = 1; rec <= AnIterTraceNumRec; rec++) {
980 AnIterTraceRec& lcAnIter = AnIterTrace[rec];
981 int toIter = lcAnIter.AnIterTraceIter;
982 double toTime = lcAnIter.AnIterTraceTime;
983 int dlIter = toIter - fmIter;
984 if (rec < AnIterTraceNumRec && dlIter != AnIterTraceIterDl)
985 printf("STRANGE: %d = dlIter != AnIterTraceIterDl = %d\n", dlIter,
986 AnIterTraceIterDl);
987 double dlTime = toTime - fmTime;
988 int iterSpeed = 0;
989 if (dlTime > 0) iterSpeed = dlIter / dlTime;
990 int lc_dual_edge_weight_mode = lcAnIter.AnIterTrace_dual_edge_weight_mode;
991 std::string str_dual_edge_weight_mode;
992 if (lc_dual_edge_weight_mode == (int)DualEdgeWeightMode::STEEPEST_EDGE)
993 str_dual_edge_weight_mode = "DSE";
994 else if (lc_dual_edge_weight_mode == (int)DualEdgeWeightMode::DEVEX)
995 str_dual_edge_weight_mode = "Dvx";
996 else if (lc_dual_edge_weight_mode == (int)DualEdgeWeightMode::DANTZIG)
997 str_dual_edge_weight_mode = "Dan";
998 else
999 str_dual_edge_weight_mode = "XXX";
1000 printf("%12d (%12d:%12d) %9.4f %12d ", dlIter, fmIter, toIter, dlTime,
1001 iterSpeed);
1002 if (report_multi) {
1003 const int pct = (100 * lcAnIter.AnIterTraceMulti);
1004 printf("| %3d ", pct);
1005 }
1006 printf("|");
1007 reportOneDensity(
1008 lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN]);
1009 reportOneDensity(
1010 lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_BTRAN_EP]);
1011 reportOneDensity(
1012 lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_PRICE_AP]);
1013 double use_row_DSE_density;
1014 if (rp_dual_steepest_edge) {
1015 if (lc_dual_edge_weight_mode ==
1016 (int)DualEdgeWeightMode::STEEPEST_EDGE) {
1017 use_row_DSE_density =
1018 lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN_DSE];
1019 } else {
1020 use_row_DSE_density = 0;
1021 }
1022 reportOneDensity(use_row_DSE_density);
1023 }
1024 printf(" | %3s ", str_dual_edge_weight_mode.c_str());
1025 if (rp_dual_steepest_edge) {
1026 double use_costly_dse;
1027 printf("| ");
1028 if (lc_dual_edge_weight_mode ==
1029 (int)DualEdgeWeightMode::STEEPEST_EDGE) {
1030 use_costly_dse = lcAnIter.AnIterTraceCostlyDse;
1031 } else {
1032 use_costly_dse = 0;
1033 }
1034 reportOneDensity(use_costly_dse);
1035 }
1036 printf("\n");
1037 fmIter = toIter;
1038 fmTime = toTime;
1039 }
1040 printf("\n");
1041 // Remove any temporary record added for the final iterations
1042 if (add_extra_record) AnIterTraceNumRec--;
1043 }
1044 }
1045 #endif
1046
iterationReport(const bool header)1047 void HighsSimplexAnalysis::iterationReport(const bool header) {
1048 if (!(iteration_report_message_level & message_level)) return;
1049 if (!header && (pivotal_row_index < 0 || entering_variable < 0)) return;
1050 reportAlgorithmPhaseIterationObjective(header,
1051 iteration_report_message_level);
1052 #ifdef HiGHSDEV
1053 reportDensity(header, iteration_report_message_level);
1054 reportIterationData(header, iteration_report_message_level);
1055 #endif
1056 HighsPrintMessage(output, message_level, iteration_report_message_level,
1057 "\n");
1058 if (!header) num_iteration_report_since_last_header++;
1059 }
1060
reportAlgorithmPhaseIterationObjective(const bool header,const int this_message_level)1061 void HighsSimplexAnalysis::reportAlgorithmPhaseIterationObjective(
1062 const bool header, const int this_message_level) {
1063 if (header) {
1064 HighsPrintMessage(output, message_level, this_message_level,
1065 " Iteration Objective ");
1066 } else {
1067 std::string algorithm;
1068 if (dualAlgorithm()) {
1069 algorithm = "Du";
1070 } else {
1071 algorithm = "Pr";
1072 }
1073 HighsPrintMessage(output, message_level, this_message_level,
1074 "%2sPh%1d %10d %20.10e", algorithm.c_str(), solve_phase,
1075 simplex_iteration_count, objective_value);
1076 }
1077 }
1078
reportInfeasibility(const bool header,const int this_message_level)1079 void HighsSimplexAnalysis::reportInfeasibility(const bool header,
1080 const int this_message_level) {
1081 if (header) {
1082 HighsPrintMessage(output, message_level, this_message_level,
1083 " Infeasibilities num(sum)");
1084 } else {
1085 if (solve_phase == 1) {
1086 HighsPrintMessage(output, message_level, this_message_level,
1087 " Ph1: %d(%g)", num_primal_infeasibilities,
1088 sum_primal_infeasibilities);
1089 } else {
1090 HighsPrintMessage(output, message_level, this_message_level,
1091 " Pr: %d(%g)", num_primal_infeasibilities,
1092 sum_primal_infeasibilities);
1093 }
1094 if (sum_dual_infeasibilities > 0) {
1095 HighsPrintMessage(output, message_level, this_message_level,
1096 "; Du: %d(%g)", num_dual_infeasibilities,
1097 sum_dual_infeasibilities);
1098 }
1099 }
1100 }
1101
1102 #ifdef HiGHSDEV
reportThreads(const bool header,const int this_message_level)1103 void HighsSimplexAnalysis::reportThreads(const bool header,
1104 const int this_message_level) {
1105 if (header) {
1106 HighsPrintMessage(output, message_level, this_message_level, " Threads");
1107 } else if (num_threads > 0) {
1108 HighsPrintMessage(output, message_level, this_message_level, " %2d|%2d|%2d",
1109 min_threads, num_threads, max_threads);
1110 } else {
1111 HighsPrintMessage(output, message_level, this_message_level, " | | ");
1112 }
1113 }
1114
reportMulti(const bool header,const int this_message_level)1115 void HighsSimplexAnalysis::reportMulti(const bool header,
1116 const int this_message_level) {
1117 if (header) {
1118 HighsPrintMessage(output, message_level, this_message_level, " Multi");
1119 } else if (average_fraction_of_possible_minor_iterations_performed >= 0) {
1120 HighsPrintMessage(
1121 output, message_level, this_message_level, " %3d%%",
1122 (int)(100 * average_fraction_of_possible_minor_iterations_performed));
1123 } else {
1124 HighsPrintMessage(output, message_level, this_message_level, " ");
1125 }
1126 }
1127
reportOneDensity(const int this_message_level,const double density)1128 void HighsSimplexAnalysis::reportOneDensity(const int this_message_level,
1129 const double density) {
1130 const int log_10_density = intLog10(density);
1131 if (log_10_density > -99) {
1132 HighsPrintMessage(output, message_level, this_message_level, " %4d",
1133 log_10_density);
1134 } else {
1135 HighsPrintMessage(output, message_level, this_message_level, " ");
1136 }
1137 }
1138
reportOneDensity(const double density)1139 void HighsSimplexAnalysis::reportOneDensity(const double density) {
1140 const int log_10_density = intLog10(density);
1141 if (log_10_density > -99) {
1142 printf(" %4d", log_10_density);
1143 } else {
1144 printf(" ");
1145 }
1146 }
1147
reportDensity(const bool header,const int this_message_level)1148 void HighsSimplexAnalysis::reportDensity(const bool header,
1149 const int this_message_level) {
1150 const bool rp_dual_steepest_edge =
1151 edge_weight_mode == DualEdgeWeightMode::STEEPEST_EDGE;
1152 if (header) {
1153 HighsPrintMessage(output, message_level, this_message_level,
1154 " C_Aq R_Ep R_Ap");
1155 if (rp_dual_steepest_edge) {
1156 HighsPrintMessage(output, message_level, this_message_level, " DSE");
1157 } else {
1158 HighsPrintMessage(output, message_level, this_message_level, " ");
1159 }
1160 } else {
1161 reportOneDensity(this_message_level, col_aq_density);
1162 reportOneDensity(this_message_level, row_ep_density);
1163 reportOneDensity(this_message_level, row_ap_density);
1164 double use_row_DSE_density;
1165 if (rp_dual_steepest_edge) {
1166 use_row_DSE_density = row_DSE_density;
1167 } else {
1168 use_row_DSE_density = 0;
1169 }
1170 reportOneDensity(this_message_level, use_row_DSE_density);
1171 }
1172 }
1173
reportInvert(const bool header,const int this_message_level)1174 void HighsSimplexAnalysis::reportInvert(const bool header,
1175 const int this_message_level) {
1176 if (header) {
1177 HighsPrintMessage(output, message_level, this_message_level, " Inv");
1178 } else {
1179 HighsPrintMessage(output, message_level, this_message_level, " %2d",
1180 invert_hint);
1181 }
1182 }
1183
1184 #ifdef HiGHSDEV
reportCondition(const bool header,const int this_message_level)1185 void HighsSimplexAnalysis::reportCondition(const bool header,
1186 const int this_message_level) {
1187 if (header) {
1188 HighsPrintMessage(output, message_level, this_message_level, " k(B)");
1189 } else {
1190 HighsPrintMessage(output, message_level, this_message_level, " %10.4g",
1191 basis_condition);
1192 }
1193 }
1194 #endif
1195
reportIterationData(const bool header,const int this_message_level)1196 void HighsSimplexAnalysis::reportIterationData(const bool header,
1197 const int this_message_level) {
1198 if (header) {
1199 HighsPrintMessage(output, message_level, this_message_level,
1200 " NumCk LvR LvC EnC DlPr "
1201 " ThDu ThPr Aa");
1202 } else {
1203 HighsPrintMessage(output, message_level, this_message_level,
1204 " %11.4g %7d %7d %7d %11.4g %11.4g %11.4g %11.4g",
1205 numerical_trouble, pivotal_row_index, leaving_variable,
1206 entering_variable, primal_delta, dual_step, primal_step,
1207 pivot_value_from_column);
1208 }
1209 }
1210
intLog10(const double v)1211 int HighsSimplexAnalysis::intLog10(const double v) {
1212 int intLog10V = -99;
1213 if (v > 0) intLog10V = log(v) / log(10.0);
1214 return intLog10V;
1215 }
1216
1217 #endif
1218
dualAlgorithm()1219 bool HighsSimplexAnalysis::dualAlgorithm() {
1220 return (simplex_strategy == SIMPLEX_STRATEGY_DUAL ||
1221 simplex_strategy == SIMPLEX_STRATEGY_DUAL_TASKS ||
1222 simplex_strategy == SIMPLEX_STRATEGY_DUAL_MULTI);
1223 }
1224