Skip to content

Commit 8c2cfeb

Browse files
committed
add simple sparse example
1 parent fc5f913 commit 8c2cfeb

File tree

1 file changed

+84
-0
lines changed

1 file changed

+84
-0
lines changed
Lines changed: 84 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,84 @@
1+
//
2+
// Copyright (c) 2022 INRIA
3+
//
4+
#include <optional>
5+
#include <random>
6+
#include <Eigen/Core>
7+
#include <proxsuite/proxqp/sparse/sparse.hpp>
8+
9+
using namespace proxsuite::proxqp;
10+
11+
int
12+
main()
13+
{
14+
15+
std::cout
16+
<< "Solve a simple example with inequality constraints using sparse ProxQP"
17+
<< std::endl;
18+
19+
// define the problem
20+
double eps_abs = 1e-9;
21+
sparse::isize dim = 3, n_eq = 0, n_in = 3, nnz = 2;
22+
23+
// random utils to generate sparse matrix
24+
std::random_device rd; // obtain a random number from hardware
25+
std::mt19937 gen(rd()); // seed the generator
26+
std::uniform_int_distribution<> int_distr(0, dim - 1); // define the range
27+
std::uniform_real_distribution<> double_distr(0, 1);
28+
29+
// cost H: first way to define sparse matrix from triplet list
30+
std::vector<Eigen::Triplet<double>>
31+
coefficients; // list of non-zeros coefficients
32+
coefficients.reserve(nnz);
33+
34+
for (sparse::isize k = 0; k < nnz; k++) {
35+
int col = int_distr(gen);
36+
int row = int_distr(gen);
37+
double val = double_distr(gen);
38+
coefficients.push_back(Eigen::Triplet<double>(col, row, val));
39+
}
40+
Eigen::SparseMatrix<double> H_spa(dim, dim);
41+
H_spa.setFromTriplets(coefficients.begin(), coefficients.end());
42+
43+
Eigen::VectorXd g = Eigen::VectorXd::Random(dim);
44+
45+
// inequality constraints C: other way to define sparse matrix from dense
46+
// matrix using sparseView
47+
Eigen::MatrixXd C = Eigen::MatrixXd(n_in, n_in);
48+
C << 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0;
49+
Eigen::SparseMatrix<double> C_spa(n_in, dim);
50+
C_spa = C.sparseView();
51+
52+
Eigen::VectorXd l = Eigen::VectorXd(dim);
53+
l << -1.0, -1.0, -1.0;
54+
55+
Eigen::VectorXd u = Eigen::VectorXd(dim);
56+
u << 1.0, 1.0, 1.0;
57+
58+
std::cout << "H:\n" << H_spa << std::endl;
59+
std::cout << "g.T:" << g.transpose() << std::endl;
60+
std::cout << "C:\n" << C_spa << std::endl;
61+
std::cout << "l.T:" << l.transpose() << std::endl;
62+
std::cout << "u.T:" << u.transpose() << std::endl;
63+
64+
// create qp object and pass some settings
65+
sparse::QP<double, long long> qp(dim, n_eq, n_in);
66+
67+
qp.settings.eps_abs = eps_abs;
68+
qp.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS;
69+
qp.settings.verbose = true;
70+
71+
// initialize qp with matrices describing the problem
72+
// note: it is also possible to use update here
73+
qp.init(H_spa, g, std::nullopt, std::nullopt, C_spa, l, u);
74+
75+
qp.solve();
76+
77+
std::cout << "primal residual: " << qp.results.info.pri_res << std::endl;
78+
std::cout << "dual residual: " << qp.results.info.dua_res << std::endl;
79+
std::cout << "total number of iteration: " << qp.results.info.iter
80+
<< std::endl;
81+
std::cout << "setup timing " << qp.results.info.setup_time << " solve time "
82+
<< qp.results.info.solve_time << std::endl;
83+
return 0;
84+
}

0 commit comments

Comments
 (0)