Skip to content

Commit 514ba38

Browse files
Yanyu000Yanyu000
committed
sad useless virtual functions
1 parent 071749d commit 514ba38

File tree

9 files changed

+1526
-133
lines changed

9 files changed

+1526
-133
lines changed

‎highs/pdlp/hipdlp/cpu_linalg.hpp‎

Lines changed: 371 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,371 @@
1+
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2+
/* */
3+
/* This file is part of the HiGHS linear optimization suite */
4+
/* */
5+
/* Available as open-source under the MIT License */
6+
/* */
7+
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
8+
/**@file pdlp/hipdlp/cpu_linalg.hpp
9+
* @brief CPU implementation of the linear algebra interface.
10+
*/
11+
#ifndef PDLP_HIPDLP_CPU_LINALG_HPP
12+
#define PDLP_HIPDLP_CPU_LINALG_HPP
13+
14+
#include "linalg.hpp"
15+
#include "Highs.h"
16+
#include <stdexcept>
17+
#include <cmath>
18+
#include <numeric>
19+
#include <algorithm>
20+
21+
namespace highs {
22+
namespace pdlp {
23+
24+
// --- Helper functions (from original linalg.cc) ---
25+
static double project_box(double x, double l, double u) {
26+
return std::max(l, std::min(x, u));
27+
}
28+
29+
static double project_non_negative(double x) { return std::max(0.0, x); }
30+
31+
// --- CPU Vector Implementation ---
32+
class CpuVector : public PdlpVector {
33+
std::vector<double> data_;
34+
public:
35+
explicit CpuVector(size_t size) : data_(size, 0.0) {}
36+
37+
void copyFromHost(const std::vector<double>& host_data) override {
38+
if (host_data.size() != data_.size()) {
39+
throw std::invalid_argument("Size mismatch in copyFromHost");
40+
}
41+
data_ = host_data;
42+
}
43+
44+
void copyToHost(std::vector<double>& host_data) const override {
45+
if (host_data.size() != data_.size()) {
46+
host_data.resize(data_.size());
47+
}
48+
host_data = data_;
49+
}
50+
void copyFrom(const PdlpVector& other) override {
51+
const auto* other_cpu = dynamic_cast<const CpuVector*>(&other);
52+
if (!other_cpu) {
53+
throw std::runtime_error("Cannot copy from non-CPU vector to CPU vector");
54+
}
55+
if (other_cpu->data_.size() != data_.size()) {
56+
throw std::runtime_error("CpuVector::copyFrom size mismatch");
57+
}
58+
data_ = other_cpu->data_;
59+
}
60+
61+
void fill(double value) override {
62+
std::fill(data_.begin(), data_.end(), value);
63+
}
64+
65+
size_t size() const override { return data_.size(); }
66+
67+
// Internal access for CpuBackend
68+
std::vector<double>& getData() { return data_; }
69+
const std::vector<double>& getData() const { return data_; }
70+
};
71+
72+
// --- CPU Sparse Matrix Implementation ---
73+
class CpuSparseMatrix : public PdlpSparseMatrix {
74+
const HighsLp* lp_; // Just holds a reference
75+
HighsSparseMatrix matrix_col_wise_;
76+
77+
public:
78+
explicit CpuSparseMatrix(const HighsLp& lp) : lp_(&lp) {
79+
// Ensure we have a column-wise copy for Ax and ATy
80+
matrix_col_wise_ = lp.a_matrix_;
81+
matrix_col_wise_.ensureColwise();
82+
}
83+
84+
size_t num_rows() const override { return matrix_col_wise_.num_row_; }
85+
size_t num_cols() const override { return matrix_col_wise_.num_col_; }
86+
87+
// Internal access
88+
const HighsSparseMatrix& getMatrix() const { return matrix_col_wise_; }
89+
const HighsLp& getLp() const { return *lp_; }
90+
};
91+
92+
// --- CPU Backend Implementation ---
93+
class CpuBackend : public LinearAlgebraBackend {
94+
private:
95+
// Helper to safely cast
96+
static CpuVector& as_cpu(PdlpVector& vec) {
97+
auto* p = dynamic_cast<CpuVector*>(&vec);
98+
if (!p) throw std::runtime_error("Invalid vector type for CpuBackend");
99+
return *p;
100+
}
101+
static const CpuVector& as_cpu(const PdlpVector& vec) {
102+
const auto* p = dynamic_cast<const CpuVector*>(&vec);
103+
if (!p) throw std::runtime_error("Invalid vector type for CpuBackend");
104+
return *p;
105+
}
106+
static const CpuSparseMatrix& as_cpu(const PdlpSparseMatrix& mat) {
107+
const auto* p = dynamic_cast<const CpuSparseMatrix*>(&mat);
108+
if (!p) throw std::runtime_error("Invalid matrix type for CpuBackend");
109+
return *p;
110+
}
111+
112+
public:
113+
std::unique_ptr<PdlpVector> createVector(size_t size) const override {
114+
return std::make_unique<CpuVector>(size);
115+
}
116+
117+
std::unique_ptr<PdlpSparseMatrix> createSparseMatrix(
118+
const HighsLp& lp) const override {
119+
return std::make_unique<CpuSparseMatrix>(lp);
120+
}
121+
122+
void Ax(PdlpVector& result, const PdlpSparseMatrix& A,
123+
const PdlpVector& x) const override {
124+
auto& res_vec = as_cpu(result).getData();
125+
const auto& mat = as_cpu(A).getMatrix();
126+
const auto& x_vec = as_cpu(x).getData();
127+
128+
std::fill(res_vec.begin(), res_vec.end(), 0.0);
129+
for (HighsInt col = 0; col < mat.num_col_; ++col) {
130+
for (HighsInt i = mat.start_[col]; i < mat.start_[col + 1]; ++i) {
131+
const HighsInt row = mat.index_[i];
132+
res_vec[row] += mat.value_[i] * x_vec[col];
133+
}
134+
}
135+
}
136+
137+
void ATy(PdlpVector& result, const PdlpSparseMatrix& A,
138+
const PdlpVector& y) const override {
139+
auto& res_vec = as_cpu(result).getData();
140+
const auto& mat = as_cpu(A).getMatrix();
141+
const auto& y_vec = as_cpu(y).getData();
142+
143+
std::fill(res_vec.begin(), res_vec.end(), 0.0);
144+
for (HighsInt col = 0; col < mat.num_col_; ++col) {
145+
for (HighsInt i = mat.start_[col]; i < mat.start_[col + 1]; ++i) {
146+
const HighsInt row = mat.index_[i];
147+
res_vec[col] += mat.value_[i] * y_vec[row];
148+
}
149+
}
150+
}
151+
152+
double dot(const PdlpVector& a, const PdlpVector& b) const override {
153+
const auto& a_vec = as_cpu(a).getData();
154+
const auto& b_vec = as_cpu(b).getData();
155+
double result = 0.0;
156+
for (size_t i = 0; i < a_vec.size(); ++i) {
157+
result += a_vec[i] * b_vec[i];
158+
}
159+
return result;
160+
}
161+
162+
double nrm2(const PdlpVector& a) const override {
163+
return std::sqrt(dot(a, a));
164+
}
165+
166+
double vector_norm(const PdlpVector& vec, double p) const override {
167+
const auto& v = as_cpu(vec).getData();
168+
if (std::isinf(p)) {
169+
double max_val = 0.0;
170+
for (double val : v) max_val = std::max(max_val, std::abs(val));
171+
return max_val;
172+
}
173+
if (p == 1.0) {
174+
double sum = 0.0;
175+
for (double val : v) sum += std::abs(val);
176+
return sum;
177+
}
178+
if (p == 2.0) {
179+
return nrm2(vec);
180+
}
181+
double sum = 0.0;
182+
for (double val : v) sum += std::pow(std::abs(val), p);
183+
return std::pow(sum, 1.0 / p);
184+
}
185+
186+
187+
void scale(PdlpVector& a, double factor) const override {
188+
auto& a_vec = as_cpu(a).getData();
189+
for (double& val : a_vec) {
190+
val *= factor;
191+
}
192+
}
193+
194+
void sub(PdlpVector& result, const PdlpVector& a,
195+
const PdlpVector& b) const override {
196+
auto& res_vec = as_cpu(result).getData();
197+
const auto& a_vec = as_cpu(a).getData();
198+
const auto& b_vec = as_cpu(b).getData();
199+
for (size_t i = 0; i < res_vec.size(); ++i) {
200+
res_vec[i] = a_vec[i] - b_vec[i];
201+
}
202+
}
203+
204+
void updateX(PdlpVector& x_new, const PdlpVector& x_old,
205+
const PdlpVector& aty, const PdlpVector& c,
206+
const PdlpVector& l, const PdlpVector& u,
207+
double primal_step) const override {
208+
auto& x_new_vec = as_cpu(x_new).getData();
209+
const auto& x_old_vec = as_cpu(x_old).getData();
210+
const auto& aty_vec = as_cpu(aty).getData();
211+
const auto& c_vec = as_cpu(c).getData();
212+
const auto& l_vec = as_cpu(l).getData();
213+
const auto& u_vec = as_cpu(u).getData();
214+
215+
for (size_t i = 0; i < x_new_vec.size(); ++i) {
216+
double gradient = c_vec[i] - aty_vec[i];
217+
double x_unproj = x_old_vec[i] - primal_step * gradient;
218+
x_new_vec[i] = project_box(x_unproj, l_vec[i], u_vec[i]);
219+
}
220+
}
221+
222+
void updateY(PdlpVector& y_new, const PdlpVector& y_old,
223+
const PdlpVector& ax, const PdlpVector& ax_next,
224+
const PdlpVector& rhs,
225+
const PdlpVector& is_equality_row_vec,
226+
double dual_step) const override {
227+
auto& y_new_vec = as_cpu(y_new).getData();
228+
const auto& y_old_vec = as_cpu(y_old).getData();
229+
const auto& ax_vec = as_cpu(ax).getData();
230+
const auto& ax_next_vec = as_cpu(ax_next).getData();
231+
const auto& rhs_vec = as_cpu(rhs).getData();
232+
const auto& is_eq_vec = as_cpu(is_equality_row_vec).getData();
233+
234+
for (size_t j = 0; j < y_new_vec.size(); j++) {
235+
double extr_ax = 2 * ax_next_vec[j] - ax_vec[j];
236+
double dual_update = y_old_vec[j] + dual_step * (rhs_vec[j] - extr_ax);
237+
bool is_equality = (is_eq_vec[j] > 0.5);
238+
y_new_vec[j] =
239+
is_equality ? dual_update : project_non_negative(dual_update);
240+
}
241+
}
242+
243+
void accumulate_weighted_sum(PdlpVector& sum_vec, const PdlpVector& vec,
244+
double weight) const override {
245+
auto& s = as_cpu(sum_vec).getData();
246+
const auto& v = as_cpu(vec).getData();
247+
for (size_t i = 0; i < s.size(); ++i) {
248+
s[i] += v[i] * weight;
249+
}
250+
}
251+
252+
void compute_average(PdlpVector& avg_vec, const PdlpVector& sum_vec,
253+
double total_weight) const override {
254+
auto& avg = as_cpu(avg_vec).getData();
255+
const auto& sum = as_cpu(sum_vec).getData();
256+
const double scale = (total_weight > 1e-10) ? 1.0 / total_weight : 1.0;
257+
for (size_t i = 0; i < avg.size(); ++i) {
258+
avg[i] = sum[i] * scale;
259+
}
260+
}
261+
262+
void computeDualSlacks(PdlpVector& dSlackPos, PdlpVector& dSlackNeg,
263+
const PdlpVector& c, const PdlpVector& aty,
264+
const PdlpVector& col_lower,
265+
const PdlpVector& col_upper,
266+
PdlpVector& dualResidual) const override {
267+
268+
auto& s_pos = as_cpu(dSlackPos).getData();
269+
auto& s_neg = as_cpu(dSlackNeg).getData();
270+
const auto& c_vec = as_cpu(c).getData();
271+
const auto& aty_vec = as_cpu(aty).getData();
272+
const auto& l_vec = as_cpu(col_lower).getData();
273+
const auto& u_vec = as_cpu(col_upper).getData();
274+
auto& res_vec = as_cpu(dualResidual).getData();
275+
276+
for (size_t i = 0; i < c_vec.size(); ++i) {
277+
res_vec[i] = c_vec[i] - aty_vec[i];
278+
279+
if (l_vec[i] > -kHighsInf) {
280+
s_pos[i] = std::max(0.0, res_vec[i]);
281+
} else {
282+
s_pos[i] = 0.0;
283+
}
284+
if (u_vec[i] < kHighsInf) {
285+
s_neg[i] = std::max(0.0, -res_vec[i]);
286+
} else {
287+
s_neg[i] = 0.0;
288+
}
289+
}
290+
}
291+
292+
double computePrimalFeasibility(
293+
const PdlpVector& ax, const PdlpVector& rhs,
294+
const PdlpVector& is_equality_row_vec,
295+
const std::vector<double>& row_scale,
296+
PdlpVector& primalResidual) const override {
297+
298+
const auto& ax_vec = as_cpu(ax).getData();
299+
const auto& rhs_vec = as_cpu(rhs).getData();
300+
const auto& is_eq_vec = as_cpu(is_equality_row_vec).getData();
301+
auto& res_vec = as_cpu(primalResidual).getData();
302+
303+
double norm_sq = 0.0;
304+
for (size_t i = 0; i < ax_vec.size(); ++i) {
305+
double residual = ax_vec[i] - rhs_vec[i];
306+
if (is_eq_vec[i] < 0.5) { // Inequality row
307+
residual = std::min(0.0, residual);
308+
}
309+
if (row_scale.size() > 0) {
310+
residual *= row_scale[i];
311+
}
312+
res_vec[i] = residual; // Store for debugging, not strictly needed
313+
norm_sq += residual * residual;
314+
}
315+
return std::sqrt(norm_sq);
316+
}
317+
318+
double computeDualFeasibility(
319+
const PdlpVector& c, const PdlpVector& aty,
320+
const PdlpVector& dSlackPos, const PdlpVector& dSlackNeg,
321+
const std::vector<double>& col_scale,
322+
PdlpVector& dualResidual) const override {
323+
324+
const auto& c_vec = as_cpu(c).getData();
325+
const auto& aty_vec = as_cpu(aty).getData();
326+
const auto& s_pos = as_cpu(dSlackPos).getData();
327+
const auto& s_neg = as_cpu(dSlackNeg).getData();
328+
auto& res_vec = as_cpu(dualResidual).getData();
329+
330+
double norm_sq = 0.0;
331+
for (size_t i = 0; i < c_vec.size(); ++i) {
332+
double residual = (c_vec[i] - aty_vec[i]) - s_pos[i] + s_neg[i];
333+
if (col_scale.size() > 0) {
334+
residual *= col_scale[i];
335+
}
336+
res_vec[i] = residual; // Store for debugging
337+
norm_sq += residual * residual;
338+
}
339+
return std::sqrt(norm_sq);
340+
}
341+
342+
void computeObjectives(
343+
const PdlpVector& c, const PdlpVector& x, const PdlpVector& rhs,
344+
const PdlpVector& y, const PdlpVector& col_lower,
345+
const PdlpVector& col_upper, const PdlpVector& dSlackPos,
346+
const PdlpVector& dSlackNeg, double offset, double& primal_obj,
347+
double& dual_obj) const override {
348+
349+
// All dot products are simple host-side loops
350+
primal_obj = dot(c, x) + offset;
351+
352+
double dual_obj_rhs = dot(rhs, y);
353+
double dual_obj_lower = dot(col_lower, dSlackPos);
354+
double dual_obj_upper = dot(col_upper, dSlackNeg);
355+
356+
dual_obj = dual_obj_rhs + dual_obj_lower - dual_obj_upper + offset;
357+
}
358+
359+
void synchronize() const override {
360+
// No-op for CPU
361+
}
362+
};
363+
364+
// Factory function implementation
365+
std::unique_ptr<LinearAlgebraBackend> createCpuBackend() {
366+
return std::make_unique<CpuBackend>();
367+
}
368+
369+
} // namespace pdlp
370+
} // namespace highs
371+
#endif

‎highs/pdlp/hipdlp/defs.hpp‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -16,8 +16,6 @@
1616

1717
#include "Highs.h"
1818

19-
enum class Device { CPU, GPU };
20-
2119
enum class ScalingMethod { NONE, RUIZ, POCK_CHAMBOLLE, L2_NORM, COMBINED };
2220

2321
enum class RestartStrategy { NO_RESTART, FIXED_RESTART, ADAPTIVE_RESTART };

0 commit comments

Comments
 (0)