Skip to content

Commit 49c9598

Browse files
jajhalljajhall
committed
Added highs/pdlp/hipdlp and linalg.cc
1 parent 7917e73 commit 49c9598

File tree

4 files changed

+258
-0
lines changed

4 files changed

+258
-0
lines changed

‎CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -734,6 +734,7 @@ if(NOT FAST_BUILD)
734734
${PROJECT_SOURCE_DIR}/highs
735735
${PROJECT_SOURCE_DIR}/highs/io
736736
${PROJECT_SOURCE_DIR}/highs/pdlp/cupdlp
737+
${PROJECT_SOURCE_DIR}/highs/pdlp/hipdlp
737738
${PROJECT_SOURCE_DIR}/highs/ipm/ipx
738739
${PROJECT_SOURCE_DIR}/highs/ipm/basiclu
739740
${PROJECT_SOURCE_DIR}/highs/lp_data

‎cmake/sources-python.cmake‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@ set(include_dirs_python
1515
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/highs/parallel>
1616
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/highs/pdlp>
1717
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/highs/pdlp/cupdlp>
18+
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/highs/pdlp/hipdlp>
1819
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/highs/presolve>
1920
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/highs/qpsolver>
2021
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/highs/simplex>
@@ -235,6 +236,7 @@ set(highs_sources_python
235236
highs/model/HighsModel.cpp
236237
highs/parallel/HighsTaskExecutor.cpp
237238
highs/pdlp/CupdlpWrapper.cpp
239+
highs/pdlp/hipdlp/linalg.cc
238240
highs/presolve/HighsPostsolveStack.cpp
239241
highs/presolve/HighsSymmetry.cpp
240242
highs/presolve/HPresolve.cpp

‎cmake/sources.cmake‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@ set(include_dirs
1515
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/highs/parallel>
1616
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/highs/pdlp>
1717
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/highs/pdlp/cupdlp>
18+
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/highs/pdlp/hipdlp>
1819
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/highs/presolve>
1920
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/highs/qpsolver>
2021
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/highs/simplex>
@@ -240,6 +241,7 @@ set(highs_sources
240241
model/HighsModel.cpp
241242
parallel/HighsTaskExecutor.cpp
242243
pdlp/CupdlpWrapper.cpp
244+
pdlp/hipdlp/linalg.cc
243245
presolve/HighsPostsolveStack.cpp
244246
presolve/HighsSymmetry.cpp
245247
presolve/HPresolve.cpp

‎highs/pdlp/hipdlp/linalg.cc‎

Lines changed: 253 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,253 @@
1+
/*
2+
* @Author: Zhou Yanyu（周妍妤） [email protected]
3+
* @Date: 2025-07-09 14:54:26
4+
* @LastEditors: Zhou Yanyu（周妍妤） [email protected]
5+
* @LastEditTime: 2025-08-05 14:47:19
6+
* @FilePath: /cupdlp-CPP/src/linalg.cpp
7+
* @Description: 这是默认设置,请设置`customMade`, 打开koroFileHeader查看配置 进行设置: https://github.com/OBKoro1/koro1FileHeader/wiki/%E9%85%8D%E7%BD%AE
8+
*/
9+
#include "linalg.hpp"
10+
#include "Highs.h"
11+
12+
namespace linalg {
13+
14+
double project_box(double x, double l, double u){
15+
if (x < l) {
16+
return l;
17+
} else if (x > u) {
18+
return u;
19+
} else {
20+
return x;
21+
}
22+
}
23+
24+
double project_non_negative(double x){
25+
return (x < 0.0) ? 0.0 : x;
26+
}
27+
28+
void Ax(const HighsLp &lp, const std::vector<double> &x, std::vector<double> & result){
29+
for (HighsInt i = 0; i < lp.num_row_; ++i) {
30+
result[i] = 0.0;
31+
}
32+
for (HighsInt col = 0; col < lp.num_col_; ++col) { // Loop over columns
33+
for (HighsInt k_ptr = lp.a_matrix_.start_[col]; k_ptr < lp.a_matrix_.start_[col + 1]; ++k_ptr) {
34+
HighsInt row = lp.a_matrix_.index_[k_ptr]; // row index
35+
double a_val = lp.a_matrix_.value_[k_ptr];
36+
result[row] += a_val * x[col]; // A[row][col] * x[col]
37+
}
38+
}
39+
}
40+
41+
void ATy(const HighsLp &lp, const std::vector<double> &y, std::vector<double> & result){
42+
for (HighsInt i = 0; i < lp.num_col_; ++i) {
43+
result[i] = 0.0;
44+
}
45+
for (HighsInt col = 0; col < lp.num_col_; ++col) { // Loop over columns
46+
for (HighsInt k_ptr = lp.a_matrix_.start_[col]; k_ptr < lp.a_matrix_.start_[col + 1]; ++k_ptr) {
47+
HighsInt row = lp.a_matrix_.index_[k_ptr]; // row index
48+
double a_val = lp.a_matrix_.value_[k_ptr]; // A[row][col]
49+
result[col] += a_val * y[row]; // A^T[col][row] * y[row]
50+
}
51+
}
52+
}
53+
54+
double nrm2(const std::vector<double>& vec) {
55+
double sum_sq = 0.0;
56+
for (double val : vec) {
57+
sum_sq += val * val;
58+
}
59+
return std::sqrt(sum_sq);
60+
}
61+
62+
// ADD THE IMPLEMENTATION FOR scale:
63+
void scale(std::vector<double>& vec, double factor) {
64+
for (size_t i = 0; i < vec.size(); ++i) {
65+
vec[i] *= factor;
66+
}
67+
}
68+
69+
void normalize(std::vector<double>& vec){
70+
double norm = nrm2(vec);
71+
if (norm > 0.0) {
72+
scale(vec, 1.0 / norm);
73+
} else {
74+
// If the vector is zero, we can choose to leave it as is or set it to a default value
75+
// Here we choose to leave it unchanged
76+
}
77+
}
78+
79+
double dot(const std::vector<double>& a, const std::vector<double>& b) {
80+
if (a.size() != b.size()) {
81+
throw std::invalid_argument("Vectors must be of the same size");
82+
}
83+
double result = 0.0;
84+
for (size_t i = 0; i < a.size(); ++i) {
85+
result += a[i] * b[i];
86+
}
87+
return result;
88+
}
89+
90+
// Computes the L2 norm of the difference between two vectors (v1 - v2).
91+
double diffTwoNorm(const std::vector<double>& v1, const std::vector<double>& v2) {
92+
double norm_sq = 0.0;
93+
if (v1.size() != v2.size()) {
94+
// Handle error: vectors must have the same dimension.
95+
return -1.0;
96+
}
97+
for (size_t i = 0; i < v1.size(); ++i) {
98+
double diff = v1[i] - v2[i];
99+
norm_sq += diff * diff;
100+
}
101+
return std::sqrt(norm_sq);
102+
}
103+
104+
double vector_norm(const std::vector<double>& vec, double p) {
105+
if (std::isinf(p)) {
106+
// Infinity norm
107+
double max_val = 0.0;
108+
for (double val : vec) {
109+
max_val = std::max(max_val, std::abs(val));
110+
}
111+
return max_val;
112+
} else if (p == 2.0) {
113+
// L2 norm (use existing optimized function)
114+
return nrm2(vec);
115+
} else if (p == 1.0) {
116+
// L1 norm
117+
double sum = 0.0;
118+
for (double val : vec) {
119+
sum += std::abs(val);
120+
}
121+
return sum;
122+
} else {
123+
// General p-norm
124+
double sum = 0.0;
125+
for (double val : vec) {
126+
sum += std::pow(std::abs(val), p);
127+
}
128+
return std::pow(sum, 1.0 / p);
129+
}
130+
}
131+
132+
// General vector norm (for raw array)
133+
double vector_norm(const double* values, size_t size, double p) {
134+
if (std::isinf(p)) {
135+
// Infinity norm
136+
double max_val = 0.0;
137+
for (size_t i = 0; i < size; ++i) {
138+
max_val = std::max(max_val, std::abs(values[i]));
139+
}
140+
return max_val;
141+
} else if (p == 2.0) {
142+
// L2 norm
143+
double sum_sq = 0.0;
144+
for (size_t i = 0; i < size; ++i) {
145+
sum_sq += values[i] * values[i];
146+
}
147+
return std::sqrt(sum_sq);
148+
} else if (p == 1.0) {
149+
// L1 norm
150+
double sum = 0.0;
151+
for (size_t i = 0; i < size; ++i) {
152+
sum += std::abs(values[i]);
153+
}
154+
return sum;
155+
} else {
156+
// General p-norm
157+
double sum = 0.0;
158+
for (size_t i = 0; i < size; ++i) {
159+
sum += std::pow(std::abs(values[i]), p);
160+
}
161+
return std::pow(sum, 1.0 / p);
162+
}
163+
}
164+
165+
double compute_cost_norm(const HighsLp& lp, double p) {
166+
return vector_norm(lp.col_cost_, p);
167+
}
168+
169+
// Compute norm of RHS vector (only for finite values)
170+
double compute_rhs_norm(const HighsLp& lp, double p) {
171+
std::vector<double> finite_rhs;
172+
finite_rhs.reserve(lp.num_row_);
173+
174+
for (HighsInt i = 0; i < lp.num_row_; ++i) {
175+
if (lp.row_lower_[i] > -kHighsInf && lp.row_lower_[i] < kHighsInf) {
176+
finite_rhs.push_back(lp.row_lower_[i]);
177+
}
178+
}
179+
180+
return finite_rhs.empty() ? 0.0 : vector_norm(finite_rhs, p);
181+
}
182+
183+
// Compute column norms of the constraint matrix
184+
std::vector<double> compute_column_norms(const HighsLp& lp, double p) {
185+
std::vector<double> col_norms(lp.num_col_, 0.0);
186+
187+
for (HighsInt col = 0; col < lp.num_col_; ++col) {
188+
HighsInt start = lp.a_matrix_.start_[col];
189+
HighsInt end = lp.a_matrix_.start_[col + 1];
190+
191+
if (start < end) {
192+
col_norms[col] = vector_norm(&lp.a_matrix_.value_[start], end - start, p);
193+
}
194+
}
195+
196+
return col_norms;
197+
}
198+
199+
// Compute row norms of the constraint matrix
200+
std::vector<double> compute_row_norms(const HighsLp& lp, double p) {
201+
std::vector<double> row_norms(lp.num_row_, 0.0);
202+
203+
if (std::isinf(p)) {
204+
// Infinity norm - find max absolute value in each row
205+
for (HighsInt col = 0; col < lp.num_col_; ++col) {
206+
for (HighsInt el = lp.a_matrix_.start_[col];
207+
el < lp.a_matrix_.start_[col + 1]; ++el) {
208+
HighsInt row = lp.a_matrix_.index_[el];
209+
double abs_val = std::abs(lp.a_matrix_.value_[el]);
210+
row_norms[row] = std::max(row_norms[row], abs_val);
211+
}
212+
}
213+
} else if (p == 2.0) {
214+
// L2 norm - sum of squares
215+
for (HighsInt col = 0; col < lp.num_col_; ++col) {
216+
for (HighsInt el = lp.a_matrix_.start_[col];
217+
el < lp.a_matrix_.start_[col + 1]; ++el) {
218+
HighsInt row = lp.a_matrix_.index_[el];
219+
double val = lp.a_matrix_.value_[el];
220+
row_norms[row] += val * val;
221+
}
222+
}
223+
for (HighsInt row = 0; row < lp.num_row_; ++row) {
224+
row_norms[row] = std::sqrt(row_norms[row]);
225+
}
226+
} else if (p == 1.0) {
227+
// L1 norm - sum of absolute values
228+
for (HighsInt col = 0; col < lp.num_col_; ++col) {
229+
for (HighsInt el = lp.a_matrix_.start_[col];
230+
el < lp.a_matrix_.start_[col + 1]; ++el) {
231+
HighsInt row = lp.a_matrix_.index_[el];
232+
row_norms[row] += std::abs(lp.a_matrix_.value_[el]);
233+
}
234+
}
235+
} else {
236+
// General p-norm
237+
for (HighsInt col = 0; col < lp.num_col_; ++col) {
238+
for (HighsInt el = lp.a_matrix_.start_[col];
239+
el < lp.a_matrix_.start_[col + 1]; ++el) {
240+
HighsInt row = lp.a_matrix_.index_[el];
241+
row_norms[row] += std::pow(std::abs(lp.a_matrix_.value_[el]), p);
242+
}
243+
}
244+
for (HighsInt row = 0; row < lp.num_row_; ++row) {
245+
if (row_norms[row] > 0.0) {
246+
row_norms[row] = std::pow(row_norms[row], 1.0 / p);
247+
}
248+
}
249+
}
250+
251+
return row_norms;
252+
}
253+
}

0 commit comments

Comments
 (0)