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