| 47 | } |
| 48 | |
| 49 | double GetBasisSolvesCheckSolution(const HighsLp& lp, |
| 50 | const vector<HighsInt>& basic_variables, |
| 51 | const vector<double>& rhs, |
| 52 | const vector<double>& solution, |
| 53 | const bool transpose = false) { |
| 54 | const double residual_tolerance = 1e-8; |
| 55 | double residual_norm = 0; |
| 56 | if (transpose) { |
| 57 | for (HighsInt k = 0; k < lp.num_row_; k++) { |
| 58 | double residual = 0; |
| 59 | HighsInt var = basic_variables[k]; |
| 60 | if (var < 0) { |
| 61 | HighsInt row = -(1 + var); |
| 62 | residual = fabs(rhs[k] - solution[row]); |
| 63 | if (residual > residual_tolerance) { |
| 64 | if (dev_run) |
| 65 | printf("Row |[B^Tx-b]_{%2" HIGHSINT_FORMAT "}| = %11.4g\n", k, |
| 66 | residual); |
| 67 | } |
| 68 | } else { |
| 69 | HighsInt col = var; |
| 70 | for (HighsInt el = lp.a_matrix_.start_[col]; |
| 71 | el < lp.a_matrix_.start_[col + 1]; el++) { |
| 72 | HighsInt row = lp.a_matrix_.index_[el]; |
| 73 | residual += lp.a_matrix_.value_[el] * solution[row]; |
| 74 | } |
| 75 | residual = fabs(rhs[k] - residual); |
| 76 | if (residual > residual_tolerance) { |
| 77 | if (dev_run) |
| 78 | printf("Col |[B^Tx-b]_{%2" HIGHSINT_FORMAT "}| = %11.4g\n", k, |
| 79 | residual); |
| 80 | } |
| 81 | } |
| 82 | residual_norm += residual; |
| 83 | } |
| 84 | } else { |
| 85 | vector<double> basis_matrix_times_solution; |
| 86 | basis_matrix_times_solution.assign(lp.num_row_, 0); |
| 87 | for (HighsInt k = 0; k < lp.num_row_; k++) { |
| 88 | HighsInt var = basic_variables[k]; |
| 89 | if (var < 0) { |
| 90 | HighsInt row = -(1 + var); |
| 91 | basis_matrix_times_solution[row] += solution[k]; |
| 92 | } else { |
| 93 | HighsInt col = var; |
| 94 | for (HighsInt el = lp.a_matrix_.start_[col]; |
| 95 | el < lp.a_matrix_.start_[col + 1]; el++) { |
| 96 | HighsInt row = lp.a_matrix_.index_[el]; |
| 97 | basis_matrix_times_solution[row] += |
| 98 | lp.a_matrix_.value_[el] * solution[k]; |
| 99 | } |
| 100 | } |
| 101 | } |
| 102 | for (HighsInt k = 0; k < lp.num_row_; k++) { |
| 103 | double residual = fabs(rhs[k] - basis_matrix_times_solution[k]); |
| 104 | if (residual > residual_tolerance) { |
| 105 | if (dev_run) |
| 106 | printf("|[B^Tx-b]_{%2" HIGHSINT_FORMAT "}| = %11.4g\n", k, residual); |