1 #include <barvinok/basis_reduction.h>
3 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
5 static void save_alpha(GBR_LP
*lp
, int first
, int n
, GBR_type
*alpha
)
9 for (i
= 0; i
< n
; ++i
)
10 GBR_lp_get_alpha(lp
, first
+i
, &alpha
[i
]);
13 /* This function implements the algorithm described in
14 * "An Implementation of the Generalized Basis Reduction Algorithm
15 * for Integer Programming" of Cook el al. to compute a reduced basis.
16 * We use \epsilon = 1/4.
18 Matrix
*Polyhedron_Reduced_Basis(Polyhedron
*P
)
20 int dim
= P
->Dimension
;
22 Matrix
*basis
= Identity(dim
);
24 GBR_type F_old
, alpha
, F_new
;
29 GBR_type
*alpha_buffer
[2];
30 GBR_type
*alpha_saved
;
36 if (P
->Dimension
== 1)
45 b_tmp
= Vector_Alloc(dim
);
47 F
= ALLOCN(GBR_type
, dim
);
48 alpha_buffer
[0] = ALLOCN(GBR_type
, dim
);
49 alpha_buffer
[1] = ALLOCN(GBR_type
, dim
);
50 alpha_saved
= alpha_buffer
[0];
52 for (i
= 0; i
< dim
; ++i
) {
54 GBR_init(alpha_buffer
[0][i
]);
55 GBR_init(alpha_buffer
[1][i
]);
68 GBR_lp_set_obj(lp
, basis
->p
[0], dim
);
71 GBR_lp_get_obj_val(lp
, &F
[0]);
75 row
= GBR_lp_next_row(lp
);
76 GBR_set(F_new
, F_saved
);
77 GBR_set(alpha
, alpha_saved
[i
]);
79 row
= GBR_lp_add_row(lp
, basis
->p
[i
], dim
);
80 GBR_lp_set_obj(lp
, basis
->p
[i
+1], dim
);
83 GBR_lp_get_obj_val(lp
, &F_new
);
85 GBR_lp_get_alpha(lp
, row
, &alpha
);
88 save_alpha(lp
, row
-i
, i
, alpha_saved
);
92 GBR_set(F
[i
+1], F_new
);
94 GBR_floor(mu
[0], alpha
);
95 GBR_ceil(mu
[1], alpha
);
97 if (value_eq(mu
[0], mu
[1]))
98 value_assign(tmp
, mu
[0]);
102 for (j
= 0; j
<= 1; ++j
) {
103 value_assign(tmp
, mu
[j
]);
104 Vector_Combine(basis
->p
[i
+1], basis
->p
[i
], b_tmp
->p
, one
, tmp
, dim
);
105 GBR_lp_set_obj(lp
, b_tmp
->p
, dim
);
106 if (GBR_lp_solve(lp
))
108 GBR_lp_get_obj_val(lp
, &mu_F
[j
]);
110 save_alpha(lp
, row
-i
, i
, alpha_buffer
[j
]);
113 if (GBR_lt(mu_F
[0], mu_F
[1]))
118 value_assign(tmp
, mu
[j
]);
119 GBR_set(F_new
, mu_F
[j
]);
120 alpha_saved
= alpha_buffer
[j
];
122 Vector_Combine(basis
->p
[i
+1], basis
->p
[i
], basis
->p
[i
+1], one
, tmp
, dim
);
124 GBR_set(F_old
, F
[i
]);
127 /* mu_F[0] = 4 * F_new; mu_F[1] = 3 * F_old */
128 GBR_set_ui(mu_F
[0], 4);
129 GBR_mul(mu_F
[0], mu_F
[0], F_new
);
130 GBR_set_ui(mu_F
[1], 3);
131 GBR_mul(mu_F
[1], mu_F
[1], F_old
);
132 if (GBR_lt(mu_F
[0], mu_F
[1])) {
133 Vector_Exchange(basis
->p
[i
], basis
->p
[i
+1], dim
);
136 GBR_set(F_saved
, F_new
);
140 GBR_set(F
[0], F_new
);
142 GBR_lp_add_row(lp
, basis
->p
[i
], dim
);
158 for (i
= 0; i
< dim
; ++i
) {
160 GBR_clear(alpha_buffer
[0][i
]);
161 GBR_clear(alpha_buffer
[1][i
]);
164 free(alpha_buffer
[0]);
165 free(alpha_buffer
[1]);