1 #include "problem_utils.h"
2 #include "scs.h"
3 #include "scs_matrix.h"
4 #include <time.h> /* to seed random */
5
6 /*
7 create data for problem:
8
9 minimize c'*x
10 subject to Ax <=_K b
11
12 where K is a product of zero, linear, and second-order cones. A is a sparse
13 matrix in
14 CSC format. A is factor * n by n with about sqrt(n) nonzeros per column.
15
16 Construct data in such a way that the problem is primal and dual
17 feasible and thus bounded.
18 */
19
main(int argc,char ** argv)20 int main(int argc, char **argv) {
21 scs_int n, m, col_nnz, nnz, i, q_total, q_num_rows, max_q;
22 ScsCone *k;
23 ScsData *d;
24 ScsSettings *stgs;
25 ScsSolution *sol, *opt_sol;
26 ScsInfo info = {0};
27 scs_float p_f, p_l;
28 scs_int factor = 4;
29 int seed = 0;
30
31 /* default parameters */
32 p_f = 0.1;
33 p_l = 0.3;
34 seed = time(SCS_NULL);
35
36 switch (argc) {
37 case 5:
38 seed = atoi(argv[4]);
39 /* no break */
40 case 4:
41 p_f = atof(argv[2]);
42 p_l = atof(argv[3]);
43 /* no break */
44 case 2:
45 n = atoi(argv[1]);
46 break;
47 default:
48 scs_printf("usage:\t%s n p_f p_l s\n"
49 "\tcreates an SOCP with n variables where p_f fraction of "
50 "rows correspond\n"
51 "\tto equality constraints, p_l fraction of rows correspond "
52 "to LP constraints,\n"
53 "\tand the remaining percentage of rows are involved in "
54 "second-order\n"
55 "\tcone constraints. the random number generator is seeded "
56 "with s.\n"
57 "\tnote that p_f + p_l should be less than or equal to 1, "
58 "and that\n"
59 "\tp_f should be less than .33, since that corresponds to "
60 "as many equality\n"
61 "\tconstraints as variables.\n",
62 argv[0]);
63 scs_printf("\nusage:\t%s n p_f p_l\n"
64 "\tdefaults the seed to the system time\n",
65 argv[0]);
66 scs_printf("\nusage:\t%s n\n"
67 "\tdefaults to using p_f = 0.1 and p_l = 0.3\n",
68 argv[0]);
69 return 0;
70 }
71 scs_printf("seed : %i\n", seed);
72
73 k = (ScsCone *)scs_calloc(1, sizeof(ScsCone));
74 d = (ScsData *)scs_calloc(1, sizeof(ScsData));
75 stgs = (ScsSettings *)scs_calloc(1, sizeof(ScsSettings));
76 sol = (ScsSolution *)scs_calloc(1, sizeof(ScsSolution));
77 opt_sol = (ScsSolution *)scs_calloc(1, sizeof(ScsSolution));
78
79 m = factor * n;
80 col_nnz = (int)ceil(sqrt(n));
81 nnz = n * col_nnz;
82
83 max_q = (scs_int)ceil(factor * n / log(factor * n));
84
85 if (p_f + p_l > 1.0) {
86 printf("error: p_f + p_l > 1.0!\n");
87 return 1;
88 }
89
90 k->z = (scs_int)floor(factor * n * p_f);
91 k->l = (scs_int)floor(factor * n * p_l);
92
93 k->qsize = 0;
94 q_num_rows = factor * n - k->z - k->l;
95 k->q = (scs_int *)scs_malloc(q_num_rows * sizeof(scs_int));
96
97 while (q_num_rows > max_q) {
98 int size;
99 size = (rand() % max_q) + 1;
100 k->q[k->qsize] = size;
101 k->qsize++;
102 q_num_rows -= size;
103 }
104 if (q_num_rows > 0) {
105 k->q[k->qsize] = q_num_rows;
106 k->qsize++;
107 }
108
109 q_total = 0;
110 for (i = 0; i < k->qsize; i++) {
111 q_total += k->q[i];
112 }
113
114 k->s = SCS_NULL;
115 k->ssize = 0;
116 k->ep = 0;
117 k->ed = 0;
118 k->bsize = 0;
119 k->bu = SCS_NULL;
120 k->bl = SCS_NULL;
121
122 scs_printf("\nA is %ld by %ld, with %ld nonzeros per column.\n", (long)m,
123 (long)n, (long)col_nnz);
124 scs_printf("A has %ld nonzeros (%f%% dense).\n", (long)nnz,
125 100 * (scs_float)col_nnz / m);
126 scs_printf("Nonzeros of A take %f GB of storage.\n",
127 ((scs_float)nnz * sizeof(scs_float)) / POWF(2, 30));
128 scs_printf("Row idxs of A take %f GB of storage.\n",
129 ((scs_float)nnz * sizeof(scs_int)) / POWF(2, 30));
130 scs_printf("Col ptrs of A take %f GB of storage.\n\n",
131 ((scs_float)n * sizeof(scs_int)) / POWF(2, 30));
132
133 printf("ScsCone information:\n");
134 printf("Zero cone rows: %ld\n", (long)k->z);
135 printf("LP cone rows: %ld\n", (long)k->l);
136 printf("Number of second-order cones: %ld, covering %ld rows, with sizes\n[",
137 (long)k->qsize, (long)q_total);
138 for (i = 0; i < k->qsize; i++) {
139 printf("%ld, ", (long)k->q[i]);
140 }
141 printf("]\n");
142 printf("Number of rows covered is %ld out of %ld.\n\n",
143 (long)(q_total + k->z + k->l), (long)m);
144
145 /* set up SCS structures */
146 d->m = m;
147 d->n = n;
148 gen_random_prob_data(nnz, col_nnz, d, k, opt_sol, seed);
149 SCS(set_default_settings)(stgs);
150
151 /* stgs->write_data_filename = "random_socp_prob"; */
152
153 scs_printf("true pri opt = %4f\n", SCS(dot)(d->c, opt_sol->x, d->n));
154 scs_printf("true dua opt = %4f\n", -SCS(dot)(d->b, opt_sol->y, d->m));
155 /* solve! */
156 scs(d, k, stgs, sol, &info);
157 scs_printf("true pri opt = %4f\n", SCS(dot)(d->c, opt_sol->x, d->n));
158 scs_printf("true dua opt = %4f\n", -SCS(dot)(d->b, opt_sol->y, d->m));
159
160 if (sol->x) {
161 scs_printf("scs pri obj= %4f\n", SCS(dot)(d->c, sol->x, d->n));
162 }
163 if (sol->y) {
164 scs_printf("scs dua obj = %4f\n", -SCS(dot)(d->b, sol->y, d->m));
165 }
166
167 SCS(free_data)(d, k, stgs);
168 SCS(free_sol)(sol);
169 SCS(free_sol)(opt_sol);
170
171 return 0;
172 }
173