Line data Source code
1 : #include <math.h>
2 : #include <stdio.h>
3 : #include <stdlib.h>
4 : #include <string.h>
5 :
6 : typedef struct simplex {
7 : int m;
8 : int n;
9 : int *var;
10 : double **a;
11 : double *b;
12 : double *x;
13 : double *c;
14 : double y;
15 : } simplex_t;
16 :
17 : typedef struct node {
18 : int m;
19 : int n;
20 : int k;
21 : int h;
22 : double xh;
23 : double ak;
24 : double bk;
25 : double *min;
26 : double *max;
27 : double **a;
28 : double *b;
29 : double *x;
30 : double *c;
31 : double z;
32 : } node_t;
33 :
34 : typedef struct link {
35 : node_t *data;
36 : struct link *next;
37 : } link_t;
38 :
39 1437803 : void free_node(node_t **p) {
40 1437803 : if (*p == NULL)
41 18960 : return;
42 1418843 : if ((*p)->a != NULL) {
43 14403332 : for (int i = 0; i < (*p)->m + 1; i++) {
44 13741808 : free((*p)->a[i]);
45 13741808 : (*p)->a[i] = NULL;
46 : }
47 661524 : free((*p)->a);
48 661524 : (*p)->a = NULL;
49 : }
50 1418843 : if ((*p)->b != NULL) {
51 661524 : free((*p)->b);
52 661524 : (*p)->b = NULL;
53 : }
54 1418843 : if ((*p)->c != NULL) {
55 661524 : free((*p)->c);
56 661524 : (*p)->c = NULL;
57 : }
58 1418843 : if ((*p)->x != NULL) {
59 661524 : free((*p)->x);
60 661524 : (*p)->x = NULL;
61 : }
62 1418843 : if ((*p)->min != NULL) {
63 1418843 : free((*p)->min);
64 1418843 : (*p)->min = NULL;
65 : }
66 1418843 : if ((*p)->max != NULL) {
67 1418843 : free((*p)->max);
68 1418843 : (*p)->max = NULL;
69 : }
70 1418843 : free(*p);
71 1418843 : *p = NULL;
72 : }
73 :
74 : double intopt(int m, int n, double **a, double *b, double *c, double *x);
75 :
76 : double simplex(int m, int n, double **a, double *b, double *c, double *x,
77 : double y);
78 :
79 : double xsimplex(int m, int n, double **a, double *b, double *c, double *x,
80 : double y, int *var, int h);
81 :
82 : void pivot(simplex_t *s, int row, int col);
83 :
84 : void prepare(simplex_t *s, int k);
85 :
86 : int initial(simplex_t *s, int m, int n, double **a, double *b, double *c,
87 : double *x, double y, int *var);
88 :
89 : int select_nonbasic(simplex_t *s);
90 :
91 : int init(simplex_t *s, int m, int n, double **a, double *b, double *c,
92 : double *x, double y, int *var);
93 :
94 : const double epsilon = 1e-9;
95 :
96 : double **make_matrix(int m, int n);
97 :
98 : node_t *initial_node(int m, int n, double **a, double *b, double *c);
99 : node_t *extend(node_t *p, int m, int n, double **a, double *b, double *c, int k,
100 : double ak, double bk);
101 : int is_integer(double *xp);
102 : int integer(node_t *p);
103 : void bound(node_t *p, link_t **h, double *zp, double *x);
104 : int branch(node_t *q, double z);
105 : void succ(node_t *p, link_t **h, int m, int n, double **a, double *b, double *c,
106 : int k, double ak, double bk, double *zp, double *x);
107 :
108 : // Lab 4
109 :
110 2612 : node_t *initial_node(int m, int n, double **a, double *b, double *c) {
111 :
112 2612 : node_t *p = calloc(1, sizeof(node_t));
113 2612 : p->a = make_matrix(m + 1, n + 1);
114 2612 : p->b = calloc(m + 1, sizeof(double));
115 2612 : p->c = calloc(n + 1, sizeof(double));
116 2612 : p->x = calloc(n + 1, sizeof(double));
117 2612 : p->min = calloc(n, sizeof(double));
118 2612 : p->max = calloc(n, sizeof(double));
119 2612 : p->m = m;
120 2612 : p->n = n;
121 :
122 : // Copy a, b, and c to p
123 : // so iterate through each element of the matrix and copy it over
124 : // hopefully the compiler can unroll this loop or SIMD it.
125 : // check assembly to optimize it
126 18388 : for (size_t i = 0; i < m; i++) {
127 138880 : for (size_t j = 0; j < n; j++) {
128 123104 : p->a[i][j] = a[i][j];
129 : }
130 15776 : p->b[i] = b[i];
131 : }
132 18388 : for (size_t k = 0; k < n; k++) {
133 15776 : p->c[k] = c[k];
134 : }
135 18388 : for (size_t i = 0; i < n; i++) {
136 15776 : p->min[i] = -INFINITY;
137 15776 : p->max[i] = INFINITY;
138 : }
139 2612 : return p;
140 : }
141 :
142 1513618 : node_t *extend(node_t *p, int m, int n, double **a, double *b, double *c, int k,
143 : double ak, double bk) {
144 :
145 1513618 : node_t *q = malloc(sizeof(*q));
146 :
147 1513618 : q->k = k;
148 1513618 : q->ak = ak;
149 1513618 : q->bk = bk;
150 :
151 1513618 : if (ak > 0 && isfinite(p->max[k])) {
152 267510 : q->m = p->m;
153 1246108 : } else if (ak < 0 && p->min[k] > 0) {
154 273296 : q->m = p->m;
155 : } else {
156 972812 : q->m = p->m + 1;
157 : }
158 :
159 1513618 : q->n = p->n;
160 1513618 : q->h = -1;
161 1513618 : q->a = make_matrix(q->m + 1, q->n + 1);
162 1513618 : q->b = calloc(q->m + 1, sizeof(double));
163 1513618 : q->c = calloc(q->n + 1, sizeof(double));
164 1513618 : q->x = calloc(q->n + 1, sizeof(double));
165 1513618 : q->min = calloc(n, sizeof(double));
166 1513618 : q->max = calloc(n, sizeof(double));
167 :
168 : // copy values to q
169 18538964 : for (size_t i = 0; i < n; i++) {
170 17025346 : q->min[i] = p->min[i];
171 17025346 : q->max[i] = p->max[i];
172 17025346 : q->c[i] = c[i];
173 : }
174 18538964 : for (size_t i = 0; i < m; i++) {
175 216730748 : for (size_t j = 0; j < n; j++)
176 199705402 : q->a[i][j] = a[i][j];
177 :
178 17025346 : q->b[i] = b[i];
179 : }
180 :
181 1513618 : if (ak > 0) {
182 756809 : if (q->max[k] == INFINITY || bk < q->max[k])
183 756809 : q->max[k] = bk;
184 :
185 756809 : } else if (q->min[k] == -INFINITY || -bk > q->min[k])
186 659484 : q->min[k] = -bk;
187 :
188 18538964 : for (size_t i = m, j = 0; j < n; j++) {
189 17025346 : if (q->min[j] > -INFINITY) {
190 5672017 : q->a[i][j] = -1;
191 5672017 : q->b[i] = -q->min[j];
192 5672017 : i++;
193 : }
194 17025346 : if (q->max[j] < INFINITY) {
195 7730721 : q->a[i][j] = 1;
196 7730721 : q->b[i] = q->max[j];
197 7730721 : i++;
198 : }
199 : }
200 :
201 1513618 : return q;
202 : }
203 :
204 11240364 : int is_integer(double *xp) {
205 11240364 : double x = *xp;
206 11240364 : double r = lround(x); // round or lround
207 11240364 : if (fabs(r - x) < epsilon) {
208 9345797 : *xp = r;
209 9345797 : return 1;
210 : }
211 1894567 : return 0;
212 : }
213 :
214 1059894 : int integer(node_t *p) {
215 6184663 : for (size_t i = 0; i < p->n; i++) {
216 6164631 : if (!is_integer(&(p->x[i]))) {
217 1039862 : return 0;
218 : }
219 : }
220 20032 : return 1;
221 : }
222 :
223 :
224 18960 : void bound(node_t *p, link_t **h, double *zp, double *x)
225 : {
226 18960 : if (p == NULL || h == NULL)
227 0 : return;
228 :
229 18960 : if (p->z > *zp)
230 : {
231 13245 : *zp = p->z;
232 13245 : memcpy(x, p->x, sizeof(double) * p->n);
233 :
234 13245 : if ((*h) == NULL)
235 : {
236 838 : return;
237 : }
238 : link_t *q, *prev, *next;
239 12407 : q = *h;
240 15735 : while (q->data->z < p->z)
241 : {
242 3601 : q = q->next;
243 3601 : if (q == NULL)
244 : {
245 273 : return;
246 : }
247 3328 : if (q->data == NULL)
248 : {
249 0 : return;
250 : }
251 : }
252 :
253 12134 : prev = q;
254 12134 : q = q->next;
255 299940 : while (q != NULL)
256 : {
257 287806 : next = q->next;
258 287806 : if (q->data->z < p->z)
259 : {
260 511 : prev->next = q->next;
261 511 : free_node(&q->data);
262 511 : free(q);
263 : }
264 : else
265 : {
266 287295 : prev = q;
267 : }
268 287806 : q = next;
269 : }
270 : }
271 : }
272 :
273 1039862 : int branch(node_t *q, double z) {
274 : double min, max;
275 :
276 1039862 : if (q->z < z) {
277 185157 : return 0;
278 : }
279 :
280 5075733 : for (size_t h = 0; h < q->n; h++) {
281 5075733 : if (!is_integer(&q->x[h])) {
282 854705 : if (q->min[h] == -INFINITY) {
283 581312 : min = 0;
284 : } else {
285 273393 : min = q->min[h];
286 : }
287 854705 : max = q->max[h];
288 :
289 854705 : if (floorf(q->x[h] < min || ceilf(q->x[h]) > max)) {
290 0 : continue;
291 : }
292 854705 : q->h = h;
293 854705 : q->xh = q->x[h];
294 :
295 : // delete each of a,b,c,x
296 854705 : if (q->a != NULL) { // delete matrix
297 19072960 : for (size_t i = 0; i < q->m + 1; i++) {
298 18218255 : free(q->a[i]);
299 18218255 : q->a[i] = NULL;
300 : }
301 854705 : free(q->a);
302 854705 : q->a = NULL;
303 : }
304 854705 : if (q->b != NULL) { // etc..
305 854705 : free(q->b);
306 854705 : q->b = NULL;
307 : }
308 854705 : if (q->c != NULL) {
309 854705 : free(q->c);
310 854705 : q->c = NULL;
311 : }
312 854705 : if (q->x != NULL) {
313 854705 : free(q->x);
314 854705 : q->x = NULL;
315 : }
316 854705 : return 1;
317 : }
318 : }
319 0 : return 0;
320 : }
321 :
322 1513618 : void succ(node_t *p, link_t **h, int m, int n, double **a, double *b, double *c,
323 : int k, double ak, double bk, double *zp, double *x) {
324 :
325 1513618 : node_t *q = extend(p, m, n, a, b, c, k, ak, bk);
326 :
327 1513618 : if (q == NULL) {
328 852629 : return;
329 : }
330 :
331 1513618 : q->z = simplex(q->m, q->n, q->a, q->b, q->c, q->x, 0);
332 :
333 1513617 : if (isfinite(q->z)) {
334 1056746 : if (integer(q)) {
335 18960 : bound(q, h, zp, x);
336 18960 : free_node(&q);
337 1037786 : } else if (branch(q, *zp)) {
338 : // add q to the linked list h
339 852629 : link_t *node = malloc(sizeof(link_t)); // free this
340 852629 : node->data = q;
341 852629 : node->next = NULL;
342 :
343 : // insert the new node
344 852629 : if (h != NULL) {
345 852629 : node->next = *h;
346 852629 : *h = node;
347 : } else {
348 0 : *h = node;
349 : }
350 852629 : return;
351 : }
352 : }
353 660988 : free_node(&q);
354 : }
355 :
356 : // Prev labs
357 0 : void print(simplex_t *s) {
358 0 : printf("%14s = %10d\n%14s = %10d\n", "m", s->m, "n", s->n);
359 0 : printf("%14s = ", "max z");
360 0 : for (size_t i = 0; i < s->m; i++) {
361 0 : printf("%10.3lf*x_%zu", s->c[i], i);
362 0 : if (i != s->m - 1) {
363 0 : printf(" + ");
364 : }
365 : }
366 0 : printf("\n");
367 0 : for (size_t i = 0; i < s->m; i++) {
368 0 : for (size_t j = 0; j < s->n; ++j) {
369 0 : printf("%10.3lf*x_%zu", s->a[i][j], j);
370 0 : if (j != s->n - 1) {
371 0 : printf(" + ");
372 : } else {
373 0 : printf(" \u2264 %10.3lf", s->b[i]);
374 : }
375 : }
376 0 : printf("\n");
377 : }
378 0 : }
379 :
380 1516230 : double **make_matrix(int m, int n) {
381 : double **a;
382 :
383 1516230 : a = calloc(m, sizeof(double *));
384 33476320 : for (int i = 0; i < m; i += 1) {
385 31960090 : a[i] = calloc(n, sizeof(double));
386 : }
387 1516230 : return a;
388 : }
389 :
390 0 : void fill_matrix(double **a, int m, int n) {
391 : int i, j;
392 0 : for (i = 0; i < m; i++)
393 0 : for (j = 0; j < n; j++)
394 0 : scanf("%lf", &a[i][j]);
395 0 : }
396 :
397 0 : void free_matrix(double **a, int m) {
398 : int i;
399 0 : for (i = 0; i < m; i++)
400 0 : free(a[i]);
401 0 : free(a);
402 0 : }
403 :
404 1514307 : void prepare(simplex_t *s, int k) {
405 1514307 : int m = s->m;
406 1514307 : int n = s->n;
407 : int i;
408 :
409 31948057 : for (i = m + n; i > n; i--)
410 30433750 : s->var[i] = s->var[i - 1];
411 :
412 1514307 : s->var[n] = m + n;
413 1514307 : n += 1;
414 :
415 31948057 : for (i = 0; i < m; i++)
416 30433750 : s->a[i][n - 1] = -1;
417 :
418 1514307 : s->x = calloc(m + n, sizeof(double)); // doesn't free the old?
419 1514307 : s->c = calloc(n, sizeof(double));
420 1514307 : s->c[n - 1] = -1;
421 1514307 : s->n = n;
422 1514307 : pivot(s, k, n - 1);
423 1514307 : }
424 :
425 1516230 : double simplex(int m, int n, double **a, double *b, double *c, double *x,
426 : double y) {
427 1516230 : return xsimplex(m, n, a, b, c, x, y, NULL, 0);
428 : }
429 :
430 3030537 : double xsimplex(int m, int n, double **a, double *b, double *c, double *x,
431 : double y, int *var, int h) {
432 : simplex_t s;
433 : int i, row, col;
434 :
435 3030537 : if (!initial(&s, m, n, a, b, c, x, y, var)) {
436 456871 : free(s.var);
437 456871 : s.var = NULL;
438 456871 : return NAN;
439 : }
440 :
441 47576106 : while ((col = select_nonbasic(&s)) >= 0) {
442 45002442 : row = -1;
443 1026991910 : for (i = 0; i < m; i++) {
444 981989468 : if (a[i][col] > epsilon &&
445 506521762 : (row < 0 || b[i] / a[i][col] < b[row] / a[row][col])) {
446 139501126 : row = i;
447 : }
448 : }
449 45002442 : if (row < 0) {
450 0 : free(s.var);
451 0 : s.var = NULL;
452 0 : return NAN;
453 : }
454 45002442 : pivot(&s, row, col);
455 : }
456 :
457 2573664 : if (h == 0) {
458 13117379 : for (i = 0; i < n; i++) {
459 12058021 : if (s.var[i] < n) {
460 2493444 : x[s.var[i]] = 0;
461 : }
462 : }
463 22314424 : for (i = 0; i < m; i++) {
464 21255066 : if (s.var[n + i] < n) {
465 9564577 : x[s.var[n + i]] = s.b[i];
466 : }
467 : }
468 1059358 : free(s.var);
469 1059358 : s.var = NULL;
470 : } else {
471 20062038 : for (i = 0; i < n; i++) {
472 18547732 : x[i] = 0;
473 : }
474 31948030 : for (i = n; i < n + m; i++) {
475 30433724 : x[i] = s.b[i - n];
476 : }
477 : }
478 2573664 : return s.y;
479 : }
480 :
481 47265334 : void pivot(simplex_t *s, int row, int col) {
482 47265334 : double **a = s->a;
483 47265334 : double *b = s->b;
484 47265334 : double *c = s->c;
485 47265334 : int m = s->m;
486 47265334 : int n = s->n;
487 : int i, j, t;
488 :
489 47265334 : t = s->var[col];
490 47265334 : s->var[col] = s->var[n + row];
491 47265334 : s->var[n + row] = t;
492 47265334 : s->y = s->y + c[col] * b[row] / a[row][col];
493 :
494 652710353 : for (i = 0; i < n; ++i) {
495 605445019 : if (i != col) {
496 558179685 : c[i] = c[i] - c[col] * a[row][i] / a[row][col];
497 : }
498 : }
499 47265334 : c[col] = -c[col] / a[row][col];
500 :
501 1075911068 : for (i = 0; i < m; ++i) {
502 1028645734 : if (i != row) {
503 981380400 : b[i] = b[i] - a[i][col] * b[row] / a[row][col];
504 : }
505 : }
506 :
507 1075911066 : for (i = 0; i < m; ++i) {
508 1028645732 : if (i != row) {
509 13885327997 : for (j = 0; j < n; ++j) {
510 12903947598 : if (j != col) {
511 11922567199 : a[i][j] = a[i][j] - a[i][col] * a[row][j] / a[row][col];
512 : }
513 : }
514 : }
515 : }
516 :
517 1075911042 : for (i = 0; i < m; ++i) {
518 1028645708 : if (i != row) {
519 981380375 : a[i][col] = -a[i][col] / a[row][col];
520 : }
521 : }
522 :
523 652710338 : for (i = 0; i < n; ++i) {
524 605445004 : if (i != col) {
525 558179671 : a[row][i] = a[row][i] / a[row][col];
526 : }
527 : }
528 :
529 47265334 : b[row] = b[row] / a[row][col];
530 47265334 : a[row][col] = 1 / a[row][col];
531 :
532 47265334 : return;
533 : }
534 :
535 3030537 : int initial(simplex_t *s, int m, int n, double **a, double *b, double *c,
536 : double *x, double y, int *var) {
537 3030537 : int i, j, k = 0;
538 : double w;
539 3030537 : k = init(s, m, n, a, b, c, x, y, var);
540 :
541 3030537 : if (b[k] >= 0) {
542 1516230 : return 1; // lab 2
543 : }
544 1514307 : prepare(s, k);
545 1514307 : n = s->n;
546 1514307 : s->y = xsimplex(m, n, s->a, s->b, s->c, s->x, 0, s->var, 1);
547 :
548 31133833 : for (i = 0; i < m + n; i++) {
549 31133833 : if (s->var[i] == m + n - 1) {
550 1514306 : if (fabs(s->x[i]) > epsilon) {
551 456871 : free(s->x);
552 456871 : free(s->c);
553 456871 : s->x = NULL;
554 456871 : s->c = NULL;
555 456871 : return 0;
556 : } else
557 1057435 : break;
558 : }
559 : }
560 :
561 1057435 : if (i >= n) {
562 10437855 : for (j = 0, k = 0; k < n; k++) {
563 9689270 : if (fabs(s->a[i - n][k]) > fabs(s->a[i - n][j]))
564 1667692 : j = k;
565 : }
566 748585 : pivot(s, i - n, j);
567 748585 : i = j;
568 : }
569 :
570 1057435 : if (i < n - 1) {
571 994507 : k = s->var[i];
572 994507 : s->var[i] = s->var[n - 1];
573 994507 : s->var[n - 1] = k;
574 21092283 : for (k = 0; k < m; k++) {
575 20097776 : w = s->a[k][n - 1];
576 20097776 : s->a[k][n - 1] = s->a[k][i];
577 20097776 : s->a[k][i] = s->a[k][i];
578 20097776 : s->a[k][i] = w;
579 : }
580 : }
581 1057435 : free(s->c);
582 1057435 : s->c = c;
583 1057435 : s->y = y;
584 22302391 : for (k = n - 1; k < n + m - 1; k++)
585 21244956 : s->var[k] = s->var[k + 1];
586 1057435 : n = s->n = s->n - 1;
587 1057435 : double *t = calloc(n, sizeof(double *));
588 :
589 13107774 : for (k = 0; k < n; k++) {
590 138298244 : for (j = 0; j < n; j++) {
591 129652032 : if (k == s->var[j]) {
592 3404127 : t[j] = t[j] + s->c[k];
593 3404127 : goto next_k;
594 : }
595 : }
596 90619678 : for (j = 0; j < m; j++) {
597 90619678 : if (s->var[n + j] == k)
598 8646212 : break;
599 : }
600 8646212 : s->y = s->y + s->c[k] * s->b[j];
601 111733296 : for (i = 0; i < n; i++)
602 103087084 : t[i] = t[i] - s->c[k] * s->a[j][i];
603 :
604 12050339 : next_k:;
605 : }
606 13107774 : for (i = 0; i < n; i++)
607 12050339 : s->c[i] = t[i];
608 1057435 : free(t);
609 1057435 : free(s->x);
610 1057435 : t = NULL;
611 1057435 : s->x = NULL;
612 1057435 : return 1;
613 : }
614 :
615 47576106 : int select_nonbasic(simplex_t *s) {
616 : int i;
617 258416668 : for (i = 0; i < s->n; i++) {
618 255843004 : if (s->c[i] > epsilon) {
619 45002442 : return i;
620 : }
621 : }
622 2573664 : return -1;
623 : }
624 :
625 3030537 : int init(simplex_t *s, int m, int n, double **a, double *b, double *c,
626 : double *x, double y, int *var) {
627 : int i, k;
628 3030537 : s->m = m;
629 3030537 : s->n = n;
630 3030537 : s->a = a;
631 3030537 : s->b = b;
632 3030537 : s->c = c;
633 3030537 : s->x = x;
634 3030537 : s->y = y;
635 3030537 : s->var = var;
636 :
637 3030537 : if (s->var == NULL) {
638 1516230 : s->var = calloc(m + n + 1, sizeof(int));
639 49001212 : for (i = 0; i < m + n; i++) {
640 47484982 : s->var[i] = i;
641 : }
642 : }
643 60877610 : for (k = 0, i = 1; i < m; i++) {
644 57847073 : if (b[i] < b[k]) {
645 5631720 : k = i;
646 : }
647 : }
648 3030537 : return k;
649 : }
650 :
651 : // körs free på p och h?
652 2075 : void free_list(link_t **head) {
653 2075 : link_t *current = *head;
654 : link_t *next;
655 :
656 2075 : while (current != NULL) {
657 0 : next = current->next; // Save next pointer before we free current
658 0 : free_node(¤t->data); // Free the node_t data using your existing free_node
659 0 : free(current); // Free the link_t structure itself
660 0 : current = next;
661 : }
662 2075 : *head = NULL; // Set the head pointer to NULL after freeing
663 2075 : }
664 :
665 2612 : double intopt(int m, int n, double **a, double *b, double *c, double *x) {
666 2612 : node_t *p = initial_node(m, n, a, b, c);
667 2612 : double z = -INFINITY;
668 2612 : link_t *h = calloc(m, sizeof(link_t));
669 2612 : h->data = p;
670 2612 : h->next = NULL;
671 :
672 2612 : p->z = simplex(p->m, p->n, p->a, p->b, p->c, p->x, 0);
673 :
674 2612 : if (integer(p) || !isfinite(p->z)) {
675 536 : z = p->z;
676 536 : if (integer(p)) {
677 : // copy px to x
678 536 : memcpy(x, p->x, sizeof(double) * p->n);
679 : }
680 536 : free_node(&p);
681 536 : free(h);
682 : //free_list(&h); // this doesn't free the entire list?
683 536 : return z;
684 : }
685 :
686 2076 : branch(p, z);
687 :
688 758884 : while (h != NULL) {
689 756809 : link_t *pop = h;
690 756809 : h = pop->next;
691 756809 : node_t *p = pop->data;
692 756809 : free(pop);
693 :
694 756809 : succ(p, &h, m, n, a, b, c, p->h, 1, floorf(p->xh), &z, x);
695 756809 : succ(p, &h, m, n, a, b, c, p->h, -1, -ceilf(p->xh), &z, x);
696 756808 : free_node(&p);
697 : }
698 : //free(h);
699 2075 : free_list(&h);
700 2075 : h = NULL;
701 2075 : if (z == -INFINITY)
702 0 : return NAN;
703 : else
704 2075 : return z;
705 : }
|
