Skip to content

Commit cb5e09a

Browse files
Yanyu000Yanyu000
committed
add dSlackPosAvg
1 parent 453cc05 commit cb5e09a

File tree

8 files changed

+87
-59
lines changed

8 files changed

+87
-59
lines changed

‎check/TestPdlp.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -336,7 +336,7 @@ TEST_CASE("pdlp-restart-add-row", "[pdlp]") {
336336

337337
TEST_CASE("hi-pdlp", "[pdlp]") {
338338
std::string model =
339-
"25fv47"; //"adlittle";//"afiro";// shell// stair //25fv47 //fit2p
339+
"adlittle"; //"adlittle";//"afiro";// shell// stair //25fv47 //fit2p
340340
std::string model_file =
341341
std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
342342
Highs h;

‎highs/pdlp/CupdlpWrapper.cpp‎

Lines changed: 0 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -209,11 +209,6 @@ HighsStatus solveLpCupdlp(const HighsOptions& options, HighsTimer& timer,
209209
debugPdlpFinalSolutionLog(w->debug_pdlp_log_file_,
210210
highs_solution.col_value.data(), lp.num_col_,
211211
highs_solution.row_dual.data(), lp.num_row_);
212-
//print col_dual
213-
for (HighsInt i = 0; i < lp.num_col_; i++) {
214-
std::cout << "Column " << i
215-
<< " dual value: " << highs_solution.col_dual[i] << std::endl;
216-
}
217212
if (w->debug_pdlp_log_file_) fclose(w->debug_pdlp_log_file_);
218213
model_status = HighsModelStatus::kUnknown;
219214
highs_solution.value_valid = value_valid;

‎highs/pdlp/HiPdlpWrapper.cpp‎

Lines changed: 0 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -74,11 +74,6 @@ HighsStatus solveLpHiPdlp(const HighsOptions& options, HighsTimer& timer,
7474
debugPdlpFinalSolutionLog(pdlp.debug_pdlp_log_file_,
7575
pdlp_solution.col_value.data(), lp.num_col_,
7676
pdlp_solution.row_dual.data(), lp.num_row_);
77-
//print col_dual
78-
for (HighsInt i = 0; i < lp.num_col_; i++) {
79-
std::cout << "Column " << i
80-
<< " dual value: " << pdlp_solution.col_dual[i] << std::endl;
81-
}
8277
pdlp.closeDebugLog();
8378

8479
// Report profiling

‎highs/pdlp/cupdlp/cupdlp_solver.c‎

Lines changed: 0 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -932,12 +932,6 @@ cupdlp_retcode PDHG_Solve(const cupdlp_int* has_variables, CUPDLPwork *pdhg) {
932932
const int iter_log_between_header = 50;
933933
int iter_log_since_header = iter_log_between_header;
934934
debugPdlpIterHeaderLog(pdhg->debug_pdlp_log_file_);
935-
//print norm of pdhg->buffer2
936-
double xx = 0.0;
937-
cupdlp_twoNormSquared(pdhg, problem->data->nCols, pdhg->buffer2,
938-
&xx);
939-
cupdlp_printf("||buffer2||_2^2 = %e\n", xx);
940-
941935
for (timers->nIter = 0; timers->nIter < settings->nIterLim; ++timers->nIter) {
942936
debugPdlpIterLog(pdhg->debug_pdlp_log_file_, timers->nIter, &pdhg->debug_pdlp_data_, pdhg->stepsize->dBeta, pdhg->stepsize->dPrimalStep, pdhg->stepsize->dDualStep);
943937
PDHG_Compute_SolvingTime(pdhg);

‎highs/pdlp/hipdlp/linalg.cc‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -213,4 +213,17 @@ std::vector<double> compute_row_norms(const HighsLp& lp, double p) {
213213
return row_norms;
214214
}
215215

216+
std::vector<double> vector_subtrac(const std::vector<double>& a,
217+
const std::vector<double>& b) {
218+
if (a.size() != b.size()) {
219+
throw std::invalid_argument(
220+
"Vectors must be of the same size for subtraction.");
221+
}
222+
std::vector<double> result(a.size());
223+
for (size_t i = 0; i < a.size(); ++i) {
224+
result[i] = a[i] - b[i];
225+
}
226+
return result;
227+
}
228+
216229
} // namespace linalg

‎highs/pdlp/hipdlp/linalg.hpp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,10 @@ std::vector<double> compute_column_norms(
5252
const HighsLp& lp, double p = std::numeric_limits<double>::infinity());
5353
std::vector<double> compute_row_norms(
5454
const HighsLp& lp, double p = std::numeric_limits<double>::infinity());
55+
56+
std::vector<double> vector_subtrac(const std::vector<double>& a,
57+
const std::vector<double>& b);
58+
5559
} // namespace linalg
5660

5761
#endif // PDLP_HIPDLP_LINALG_HPP

‎highs/pdlp/hipdlp/pdhg.cc‎

Lines changed: 52 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -568,12 +568,14 @@ void PDLPSolver::solve(std::vector<double>& x, std::vector<double>& y) {
568568
// Compute residuals for current iterate
569569
bool current_converged =
570570
checkConvergence(iter, x_current_, y_current_, Ax_cache_, ATy_cache_,
571-
params_.tolerance, current_results, "[L]");
571+
params_.tolerance, current_results, "[L]",
572+
dSlackPos_, dSlackNeg_);
572573

573574
// Compute residuals for average iterate
574575
bool average_converged =
575576
checkConvergence(iter, x_avg_, y_avg_, Ax_avg, ATy_avg,
576-
params_.tolerance, average_results, "[A]");
577+
params_.tolerance, average_results, "[A]",
578+
dSlackPosAvg_, dSlackNegAvg_);
577579
hipdlpTimerStop(kHipdlpClockConvergenceCheck);
578580

579581
debugPdlpIterHeaderLog(debug_pdlp_log_file_);
@@ -598,6 +600,8 @@ void PDLPSolver::solve(std::vector<double>& x, std::vector<double>& y) {
598600
final_iter_count_ = iter;
599601
x = x_avg_;
600602
y = y_avg_;
603+
dSlackPos_ = dSlackPosAvg_;
604+
dSlackNeg_ = dSlackNegAvg_;
601605
results_ = average_results;
602606
return solveReturn(TerminationStatus::OPTIMAL);
603607
}
@@ -753,6 +757,8 @@ void PDLPSolver::initialize() {
753757
K_times_x_diff_.resize(lp_.num_row_, 0.0);
754758
dSlackPos_.resize(lp_.num_col_, 0.0);
755759
dSlackNeg_.resize(lp_.num_col_, 0.0);
760+
dSlackPosAvg_.resize(lp_.num_col_, 0.0);
761+
dSlackNegAvg_.resize(lp_.num_col_, 0.0);
756762
}
757763

758764
// Update primal weight
@@ -875,42 +881,48 @@ double PDLPSolver::computePrimalFeasibility(
875881
return linalg::vector_norm(primal_residual);
876882
}
877883

878-
void PDLPSolver::computeDualSlacks(const std::vector<double>& ATy_vector) {
884+
void PDLPSolver::computeDualSlacks(const std::vector<double>& dualResidual,
885+
std::vector<double>& dSlackPos,
886+
std::vector<double>& dSlackNeg) {
879887
// Ensure vectors are correctly sized
880-
if (dSlackPos_.size() != lp_.num_col_) dSlackPos_.resize(lp_.num_col_);
881-
if (dSlackNeg_.size() != lp_.num_col_) dSlackNeg_.resize(lp_.num_col_);
888+
if (dSlackPos.size() != lp_.num_col_) dSlackPos.resize(lp_.num_col_);
889+
if (dSlackNeg.size() != lp_.num_col_) dSlackNeg.resize(lp_.num_col_);
882890

883891
for (HighsInt i = 0; i < lp_.num_col_; ++i) {
884-
double dual_residual = lp_.col_cost_[i] - ATy_vector[i];
885-
886892
// Compute positive slack (for lower bounds)
887-
// CUPDLP: max(dual_residual, 0) * hasLower
888-
if (lp_.col_lower_[i] > -kHighsInf) {
889-
dSlackPos_[i] = std::max(0.0, dual_residual);
890-
} else {
891-
dSlackPos_[i] = 0.0;
892-
}
893-
894-
// Compute negative slack (for upper bounds)
895-
// CUPDLP: -min(dual_residual, 0) * hasUpper
896-
if (lp_.col_upper_[i] < kHighsInf) {
897-
dSlackNeg_[i] = std::max(0.0, -dual_residual);
898-
} else {
899-
dSlackNeg_[i] = 0.0;
900-
}
893+
// CUPDLP: max(dual_residual, 0) * hasLower
894+
if (lp_.col_lower_[i] > -kHighsInf) {
895+
dSlackPos[i] = std::max(0.0, dualResidual[i]);
896+
} else {
897+
dSlackPos[i] = 0.0;
898+
}
899+
900+
// Compute negative slack (for upper bounds)
901+
// CUPDLP: -min(dual_residual, 0) * hasUpper
902+
if (lp_.col_upper_[i] < kHighsInf) {
903+
dSlackNeg[i] = std::max(0.0, -dualResidual[i]);
904+
} else {
905+
dSlackNeg[i] = 0.0;
906+
}
901907
}
902908
}
903909

904910
double PDLPSolver::computeDualFeasibility(
905-
const std::vector<double>& ATy_vector) {
906-
computeDualSlacks(ATy_vector); // This updates dSlackPos_ and dSlackNeg_
911+
const std::vector<double>& ATy_vector, std::vector<double>& dSlackPos,
912+
std::vector<double>& dSlackNeg) {
913+
std::vector<double> dualResidual(lp_.num_col_, 0.0);
914+
// dualResidual = c-A'y
915+
dualResidual = linalg::vector_subtrac(lp_.col_cost_, ATy_vector);
916+
double dualResidualNorm = linalg::vector_norm(dualResidual);
917+
918+
// Call the refactored function to populate dSlackPos and dSlackNeg
919+
computeDualSlacks(dualResidual, dSlackPos, dSlackNeg);
907920

908921
std::vector<double> dual_residual(lp_.num_col_);
909922

910923
for (HighsInt i = 0; i < lp_.num_col_; ++i) {
911924
// Matching CUPDLP: c - A'y - dSlackPos + dSlackNeg
912-
dual_residual[i] =
913-
lp_.col_cost_[i] - ATy_vector[i] - dSlackPos_[i] + dSlackNeg_[i];
925+
dual_residual[i] = dualResidual[i] - dSlackPos[i] + dSlackNeg[i];
914926
}
915927

916928
// Apply scaling if needed
@@ -973,7 +985,9 @@ PDLPSolver::computeDualityGap(const std::vector<double>& x,
973985
cTx);
974986
}
975987

976-
double PDLPSolver::computeDualObjective(const std::vector<double>& y) {
988+
double PDLPSolver::computeDualObjective(
989+
const std::vector<double>& y, const std::vector<double>& dSlackPos,
990+
const std::vector<double>& dSlackNeg) {
977991
double dual_obj = lp_.offset_;
978992

979993
// Compute b'y (or rhs'y in cuPDLP notation)
@@ -984,14 +998,14 @@ double PDLPSolver::computeDualObjective(const std::vector<double>& y) {
984998
// Add contribution from lower bounds: l'*slackPos
985999
for (int i = 0; i < lp_.num_col_; ++i) {
9861000
if (lp_.col_lower_[i] > -kHighsInf) {
987-
dual_obj += lp_.col_lower_[i] * dSlackPos_[i];
1001+
dual_obj += lp_.col_lower_[i] * dSlackPos[i];
9881002
}
9891003
}
9901004

9911005
// Subtract contribution from upper bounds: u'*slackNeg
9921006
for (int i = 0; i < lp_.num_col_; ++i) {
9931007
if (lp_.col_upper_[i] < kHighsInf) {
994-
dual_obj -= lp_.col_upper_[i] * dSlackNeg_[i];
1008+
dual_obj -= lp_.col_upper_[i] * dSlackNeg[i];
9951009
}
9961010
}
9971011

@@ -1003,16 +1017,20 @@ bool PDLPSolver::checkConvergence(const int iter, const std::vector<double>& x,
10031017
const std::vector<double>& ax_vector,
10041018
const std::vector<double>& aty_vector,
10051019
double epsilon, SolverResults& results,
1006-
const char* type) {
1007-
// Compute dual slacks first
1008-
computeDualSlacks(aty_vector);
1020+
const char* type,
1021+
// Add slack vectors as non-const references
1022+
std::vector<double>& dSlackPos,
1023+
std::vector<double>& dSlackNeg) {
1024+
// computeDualSlacks is now called inside computeDualFeasibility
10091025

10101026
// Compute primal feasibility
10111027
double primal_feasibility = computePrimalFeasibility(ax_vector);
10121028
results.primal_feasibility = primal_feasibility;
10131029

10141030
// Compute dual feasibility
1015-
double dual_feasibility = computeDualFeasibility(aty_vector);
1031+
// This will populate dSlackPos and dSlackNeg
1032+
double dual_feasibility =
1033+
computeDualFeasibility(aty_vector, dSlackPos, dSlackNeg);
10161034
results.dual_feasibility = dual_feasibility;
10171035

10181036
// Compute objectives
@@ -1022,7 +1040,8 @@ bool PDLPSolver::checkConvergence(const int iter, const std::vector<double>& x,
10221040
}
10231041
results.primal_obj = primal_obj;
10241042

1025-
double dual_obj = computeDualObjective(y);
1043+
// Pass the now-populated slack vectors to computeDualObjective
1044+
double dual_obj = computeDualObjective(y, dSlackPos, dSlackNeg);
10261045
results.dual_obj = dual_obj;
10271046

10281047
// Compute duality gap
@@ -1282,8 +1301,6 @@ void PDLPSolver::unscaleSolution(std::vector<double>& x,
12821301
const std::vector<double>& col_scale = scaling_.GetColScaling();
12831302
if (!dSlackPos_.empty() && col_scale.size() == dSlackPos_.size()) {
12841303
for (size_t i = 0; i < dSlackPos_.size(); ++i) {
1285-
std::cout << "dSlackPos_ before unscale[" << i << "] = " << dSlackPos_[i]
1286-
<< std::endl;
12871304
dSlackPos_[i] *= col_scale[i];
12881305
dSlackNeg_[i] *= col_scale[i];
12891306
}

‎highs/pdlp/hipdlp/pdhg.hpp‎

Lines changed: 17 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -61,10 +61,14 @@ class PDLPSolver {
6161
void initialize();
6262
void printConstraintInfo();
6363
bool checkConvergence(const int iter, const std::vector<double>& x,
64-
const std::vector<double>& y,
65-
const std::vector<double>& ax_vector,
66-
const std::vector<double>& aty_vector, double epsilon,
67-
SolverResults& results, const char* type);
64+
const std::vector<double>& y,
65+
const std::vector<double>& ax_vector,
66+
const std::vector<double>& aty_vector,
67+
double epsilon, SolverResults& results,
68+
const char* type,
69+
// Add slack vectors as non-const references
70+
std::vector<double>& dSlackPos,
71+
std::vector<double>& dSlackNeg);
6872
void updateAverageIterates(const std::vector<double>& x,
6973
const std::vector<double>& y,
7074
const PrimalDualParams& params, int inner_iter);
@@ -96,10 +100,14 @@ class PDLPSolver {
96100
// --- Feasibility, Duality, and KKT Checks ---
97101
std::vector<double> computeLambda(const std::vector<double>& y,
98102
const std::vector<double>& ATy_vector);
99-
double computeDualObjective(const std::vector<double>& y);
103+
double computeDualObjective(const std::vector<double>& y, const std::vector<double>& dSlackPos,
104+
const std::vector<double>& dSlackNeg);
100105
double computePrimalFeasibility(const std::vector<double>& Ax_vector);
101-
void computeDualSlacks(const std::vector<double>& ATy_vector);
102-
double computeDualFeasibility(const std::vector<double>& ATy_vector);
106+
void computeDualSlacks(const std::vector<double>& dualResidual,
107+
std::vector<double>& dSlackPos,
108+
std::vector<double>& dSlackNeg);
109+
double computeDualFeasibility(const std::vector<double>& ATy_vector, std::vector<double>& dSlackPos,
110+
std::vector<double>& dSlackNeg);
103111
std::tuple<double, double, double, double, double> computeDualityGap(
104112
const std::vector<double>& x, const std::vector<double>& y,
105113
const std::vector<double>& lambda);
@@ -135,6 +143,8 @@ class PDLPSolver {
135143
int num_rejected_steps_ = 0;
136144
std::vector<double> dSlackPos_;
137145
std::vector<double> dSlackNeg_;
146+
std::vector<double> dSlackPosAvg_;
147+
std::vector<double> dSlackNegAvg_;
138148
Timer total_timer;
139149

140150
HipdlpTimer hipdlp_timer_;

0 commit comments

Comments
 (0)