1 /* This file is part of Shapes.
3 * Shapes is free software: you can redistribute it and/or modify
4 * it under the terms of the GNU General Public License as published by
5 * the Free Software Foundation, either version 3 of the License, or
8 * Shapes is distributed in the hope that it will be useful,
9 * but WITHOUT ANY WARRANTY; without even the implied warranty of
10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 * GNU General Public License for more details.
13 * You should have received a copy of the GNU General Public License
14 * along with Shapes. If not, see <http://www.gnu.org/licenses/>.
16 * Copyright 2008 Henrik Tidefelt
21 #include "basicsimplex.h"
23 #include "autoonoff.h"
26 #include <iostream> // For debugging
28 using namespace Shapes
;
30 Computation::BasicSimplex::BasicSimplex( size_t nVars
, size_t nEqns
)
31 : nVars_( nVars
), nEqns_( nEqns
),
32 nExt_( nVars
+ nEqns
),
33 a_( new double[ nEqns_
* ( nExt_
) ] ),
34 b_( new double[ nEqns_
] ),
35 c_( new double[ nExt_
] ),
36 varSet_( new size_t[ nExt_
] ),
37 nonBasicBegin_( varSet_
),
38 nonBasicEnd_( nonBasicBegin_
+ nVars_
),
39 basicBegin_( nonBasicEnd_
),
40 basicEnd_( basicBegin_
+ nEqns_
)
43 Computation::BasicSimplex::~BasicSimplex( )
52 // Call this function with preallocated memory for the result!
54 // It solves the problem of minimizing < c, *xdst > under the constraints
57 // The constraints are stored in row major order, so consecutive memory
58 // belongs to the same constraint, not the same variable.
59 // The optimum is returned in *objdst, and the return value
60 // is true if optimization finished with *objdst < objGoal.
61 // If the return value is true, then the optimum is generally not located,
62 // but only a sufficiently good point.
64 Computation::BasicSimplex::minimize( double * xdst
,
65 double * objdst
, double objGoal
,
66 const double * c
, const double * a
, const double * b
,
67 bool changeSign
) const
69 // Setup extended formulation with equality constraints and slack variables:
71 bool infeasibleStart
= false;
75 double * end
= c_
+ nVars_
;
77 const double * src
= c
;
80 for( ; dst
!= end
; ++dst
, ++src
)
87 for( ; dst
!= end
; ++dst
, ++src
)
93 for( ; dst
!= end
; ++dst
)
100 // Initialize a and b
102 // First all of a is zeroed.
103 double * end
= a_
+ nEqns_
* ( nExt_
);
105 for( ; dst
!= end
; ++dst
)
111 const double * srcb
= b
;
113 for( size_t eq
= 0; eq
< nEqns_
; ++eq
, ++srcb
, ++dstb
)
118 infeasibleStart
= true;
121 double * dst
= a_
+ eq
* ( nExt_
);
122 const double * src
= a
+ eq
* nVars_
;
123 const double * end
= src
+ nVars_
;
124 for( ; src
!= end
; ++src
, ++dst
)
133 // The basic and non-basic sets contain indices to variables.
134 // The basic set is organized so that the basic variable in equaiton i is located at ( basic + i ).
138 for( size_t * dst
= nonBasicBegin_
; dst
!= basicEnd_
; ++dst
, ++i
)
145 if( infeasibleStart
)
149 objective
= phaseTwo( xdst
, objGoal
);
153 objective
= HUGE_VAL
;
158 objective
= phaseTwo( xdst
, objGoal
);
163 *objdst
= - objective
;
169 return objective
< objGoal
;
174 Computation::BasicSimplex::maximize( double * xdst
,
175 double * objdst
, double objGoal
,
176 const double * c
, const double * a
, const double * b
) const
178 return minimize( xdst
,
185 // Function to compute a feasible tableau.
187 Computation::BasicSimplex::phaseOne( ) const
189 throw "Not implemented: Computation::BasicSimplex::phaseOne.";
192 // Function for minimizing given a feasible initial tableau.
194 Computation::BasicSimplex::phaseTwo( double * xdst
, double objGoal
) const
196 double objective
= 0;
204 std::cerr
<< "Simplex tableau:" << std::endl
;
205 for( size_t i
= 0; i
< 3 + nExt_
* 8 + 5; ++i
)
209 std::cerr
<< std::endl
;
211 for( size_t col
= 0; col
< nExt_
; ++col
)
216 for( const size_t * src
= basicBegin_
; src
!= basicEnd_
; ++src
)
235 for( const size_t * src
= nonBasicBegin_
; src
!= nonBasicEnd_
; ++src
)
254 std::cerr
<< std::endl
;
256 for( const double * src
= c_
; src
!= c_
+ nExt_
; ++src
)
259 sprintf( buf
, "%7.3f", *src
);
260 std::cerr
<< buf
<< " " ;
262 std::cerr
<< " | obj: " << objective
<< std::endl
;
263 for( size_t i
= 0; i
< 3 + nExt_
* 8 + 5; ++i
)
267 std::cerr
<< std::endl
;
269 const double * srca
= a_
;
270 for( size_t row
= 0; row
< nEqns_
; ++row
)
272 std::cerr
<< *( basicBegin_
+ row
) << " |" ;
273 const double * end
= srca
+ nExt_
;
274 for( ; srca
!= end
; ++srca
)
277 sprintf( buf
, "%7.3f", *srca
);
278 std::cerr
<< buf
<< " " ;
280 std::cerr
<< " | " << *( b_
+ row
) << std::endl
;
282 for( size_t i
= 0; i
< 3 + nExt_
* 8 + 5; ++i
)
286 std::cerr
<< std::endl
;
287 std::cerr
<< std::endl
;
290 // End of debug prints
293 if( objective
< objGoal
)
298 // Pick a pivot column. That is, the non-basic variable with the most positive coefficient
299 size_t * nonBasicCol
= 0;
302 for( size_t * src
= nonBasicBegin_
; src
!= nonBasicEnd_
; ++src
)
304 if( *( c_
+ *src
) > best
)
306 best
= *( c_
+ *src
);
311 if( nonBasicCol
== 0 )
316 size_t col
= *nonBasicCol
;
318 // Find the outgoing row
321 double thMin
= HUGE_VAL
;
322 const double * srca
= a_
+ col
;
323 const double * srcb
= b_
;
324 for( size_t i
= 0; i
< nEqns_
; ++i
, srca
+= nExt_
, ++srcb
)
328 double tmp
= *srcb
/ *srca
;
338 // There is no constraining equation.
340 // This return value is special, and does not correspond to the returned *dstx.
341 objective
= -HUGE_VAL
;
347 // Normalize the pivot
348 double * dst
= a_
+ nExt_
* row
;
349 double tmp
= 1. / *( dst
+ col
);
350 double * end
= dst
+ nExt_
;
351 for( ; dst
!= end
; ++dst
)
355 *( b_
+ row
) *= tmp
;
361 // Note that structural zeroes are kept track of using the basic/nonBasic distinction, and need not be enforced in a.
364 double * enda
= a_
+ ( nExt_
* nEqns_
);
365 const double * pivota
= a_
+ nExt_
* row
;
367 double pivotb
= *( b_
+ row
);
368 for( ; dst
!= enda
; ++dstb
)
375 double * end
= dst
+ nExt_
;
376 const double * src
= pivota
;
377 double tmp
= *( dst
+ col
);
378 for( ; dst
!= end
; ++dst
, ++src
)
382 *dstb
-= tmp
* pivotb
;
386 // "Reduce objective"
387 const double * src
= pivota
;
389 double * end
= dst
+ nExt_
;
390 double tmp
= *( dst
+ col
);
391 for( ; dst
!= end
; ++dst
, ++src
)
395 objective
-= tmp
* pivotb
;
402 *nonBasicCol
= *( basicBegin_
+ row
);
403 *( basicBegin_
+ row
) = col
;
411 double * end
= dst
+ nVars_
;
412 for( ; dst
!= end
; ++dst
)
416 const double * srcb
= b_
;
417 for( const size_t * src
= basicBegin_
; src
!= basicEnd_
; ++src
, ++srcb
)
421 *( xdst
+ *src
) = *srcb
;