項目中需要求解動態增長的線性方程組，方程個數會不斷增加，需要使用高斯消元法來求解，每個方程添加完后，已知的變量值。
代碼中分了三個類，表達式expression，方程 equation， 高斯消元求解器gauss_eliminator。
這段代碼還是有問題，如果需要最新的，請訪問我的github https://github.com/tfjiang/gauss_elimination
1 #ifndef GAUSS_ELIMINATION_H
2 #define GAUSS_ELIMINATION_H
3
4 #include <map>
5 #include <vector>
6 #include <list>
7 #include <deque>
9 #include <iostream>
10 #include <boost/unordered_map.hpp>
11 #include <boost/static_assert.hpp>
12 #include <boost/type_traits/is_same.hpp>
13 #include <zjucad/matrix/matrix.h> // a private matrix library
14
15 namespace jtf{
16 namespace algorithm{
17
18 /**
19 * @brief store an expression A * i, where A is the coeffieient
20 *
21 */
22 template <typename T>
23 class expression{
24 public:
25 expression():index(-1),coefficient(0){}
26 expression(const size_t & index_, const T & coefficient_)
27 : index(index_), coefficient(coefficient_){}
28 size_t index;
29 T coefficient;
30 /**
31 * @brief operator < is used to determin which expression ha
32 *
33 * @param other input other expression
34 * @return bool
35 */
36 bool operator < (const expression<T> & b) const{
37 return index < b.index;
38 }
39
40 /**
41 * @brief To determin whether coefficent of this expression is zeros
42 *
43 * @return bool
44 */
45 bool is_zero() const{ return fabs(coefficient) < 1e-6;}
46 };
47
48 /**
49 * @brief make expression with node_idx and value
50 *
51 * @param node_idx
52 * @param value
53 * @return expression<T>
54 */
55 template <typename T>
56 expression<T> make_expression(const size_t & node_idx, const T & value)
57 {
58 expression<T> temp(node_idx, value);
59 return temp;
60 }
61
62 /**
63 * @brief equation class
64 *
65 */
66 template <typename T>
67 class equation{
68 public:
69 typedef typename std::list<expression<T> >::const_iterator eq_const_iterator;
70 typedef typename std::list<expression<T> >::iterator eq_iterator;
71
72 equation():value_(static_cast<T>(0)){}
73
74 eq_const_iterator begin() const {return e_vec_.begin();}
75 eq_iterator begin(){return e_vec_.begin();}
76
77 eq_const_iterator end() const {return e_vec_.end();}
78 eq_iterator end(){return e_vec_.end();}
79
80 /**
81 * @brief used to standardizate the equation, sort the expressions and merge
82 * similar items and normalization to make the first item's coefficient
83 * equal 1
84 *
85 * @return int
86 */
87 int standardization(){
88 sort_equation();
89 merge_similar();
90 normalization();
91 return 0;
92 }
93
94 /**
95 * @brief update the equation with given node idx and value
96 *
97 * @param node_idx input node idx
98 * @param node_value input node value
99 * @return int
100 */
101 int update(const size_t & node_idx, const T & node_value);
102 /**
103 * @brief sort equation accronding to the expressions
104 *
105 * @return int
106 */
107 int sort_equation(){
108 e_vec_.sort();
109 return 0;
110 }
111
112 /**
113 * @brief merge similar items, WARNING. this function should be used after
114 * sorting.
115 *
116 * @return int
117 */
118 int merge_similar();
119
120 /**
121 * @brief normalize the equations: to make the first expression's coefficient
122 * equals 1
123 *
124 * @return int
125 */
126 int normalization();
127
128 /**
129 * @brief get the equation value
130 *
131 * @return const T
132 */
133 const T& get_value() const {return value_;}
134
135 /**
136 * @brief get the state of equation:
137 * if there are no expressions,
138 * if value != 0 return -1; error equation
139 * else return 0; cleared
140 * else
141 * if there is one expression return 1; calculated finial node
142 * else return 2; not finished
143 *
144 * @return int
145 */
146 int state() const;
147
148
149 /**
150 * @brief define the minus operation
151 *
152 * @param b
153 * @return equation<T>
154 */
155 equation<T> & operator -= (const equation<T> & b);
156
157 /**
158 * @brief define the output operation
159 *
160 * @param eq
161 * @return equation<T>
162 */
163 friend std::ostream& operator << (std::ostream &output,
164 const equation<T> &eq)
165 {
166 if(eq.e_vec_.empty()){
167 output << "# ------ empty expressions with value = " << eq.value() << std::endl;
168 }else{
169 output << "# ------ SUM coeff * index , val" << std::endl
170 << "# ------ ";
171 for(equation<T>::eq_const_iterator eqcit = eq.begin(); eqcit != eq.end();){
172 const expression<T> & exp = *eqcit;
173 output << exp.coefficient << "*X" << exp.index;
174 ++eqcit;
175 if(eqcit == eq.end())
176 output << " = ";
177 else
178 output << " + ";
179 }
180 output << eq.value() << std::endl;
181 }
182 return output;
183 }
184
185
186 /**
187 * @brief get the first expression idx
188 *
189 * @return size_t
190 */
191 size_t get_prim_idx() const {
192 return e_vec_.front().index;
193 }
194
195 int add_expression(const expression<T> & exp){
196 e_vec_.push_back(exp);
197 return 0;
198 }
199 T& value() {return value_;}
200 const T& value() const {return value_;}
201 std::list<expression<T> > e_vec_;
202 private:
203 T value_;
204 };
205
206 template <typename T>
207 int equation<T>::merge_similar(){
208 typedef typename std::list<expression<T> >::iterator leit;
209 leit current = e_vec_.begin();
210 leit next = current;
211 ++next;
212 while(next != e_vec_.end()){
213 if(next->index > current->index) {
214 current = next++;
215 continue;
216 }else if(next->index == current->index){
217 current->coefficient += next->coefficient;
218 e_vec_.erase(next++);
219 if(fabs(current->coefficient) < 1e-6){
220 e_vec_.erase(current++);
221 ++next;
222 }
223 }else{
224 std::cerr << "# [error] merge similar function should only be called "
225 << "after sorting." << std::endl;
226 return __LINE__;
227 }
228 }
229 return 0;
230 }
231
232 template <typename T>
233 int equation<T>::normalization()
234 {
235 const T coeff = e_vec_.front().coefficient;
236 if(fabs(coeff) < 1e-6){
237 if(e_vec_.empty()){
238 //std::cerr << "# [info] this equation is empty." << std::endl;
239 return 0;
240 }else
241 std::cerr << "# [error] this expression should be removed." << std::endl;
242 return __LINE__;
243 }
244
245 value() /= coeff;
246 for(typename std::list<expression<T> >::iterator lit = e_vec_.begin();
247 lit != e_vec_.end(); ++lit){
248 expression<T> & ep = *lit;
249 ep.coefficient /= coeff;
250 }
251 return 0;
252 }
253
254 template <typename T>
255 int equation<T>::state() const
256 {
257 if(e_vec_.empty()){
258 if(fabs(value()) < 1e-8)
259 return 0; // is cleared
260 return -1; // is conflicted
261 }else{
262 if(e_vec_.size() == 1)
263 return 1; // finial variant
264 else
265 return 2; // not ready
266 }
267 }
268
269 template <typename T>
270 equation<T> & equation<T>::operator -= (const equation<T> & b )
271 {
272 if(&b == this) {
273 e_vec_.clear();
274 value() = 0;
275 return *this;
276 }
277
278 for(typename std::list<expression<T> >::const_iterator lecit_b =
279 b.e_vec_.begin(); lecit_b != b.e_vec_.end(); ++lecit_b){
280 const expression<T> & exp = *lecit_b;
281 const size_t &node_idx = exp.index;
282 assert(fabs(exp.coefficient) > 1e-6);
283 bool is_found = false;
284
285 for(typename std::list<expression<T> >::iterator leit_a = e_vec_.begin();
286 leit_a != e_vec_.end(); ++leit_a){
287 expression<T> & exp_a = *leit_a;
288 if(exp_a.index == node_idx){
289 exp_a.coefficient -= exp.coefficient;
290 // zeros
291 if(fabs(exp_a.coefficient) < 1e-6) {
292 e_vec_.erase(leit_a);
293 is_found = true;
294 break;
295 }
296 }
297 }
298 if(!is_found){
299 e_vec_.push_back(make_expression(node_idx, -1 * exp.coefficient));
300 }
301 }
302
303 value() -= b.value();
304 sort_equation();
305 normalization();
306
307 return *this;
308 }
309
310 template <typename T>
311 int equation<T>::update(const size_t & node_idx, const T & node_value)
312 {
313 for(typename std::list<expression<T> >::iterator leit = e_vec_.begin();
314 leit != e_vec_.end();){
315 expression<T> & exp = *leit;
316 if(exp.index == node_idx){
317 value() -= exp.coefficient * node_value;
318 e_vec_.erase(leit++);
319 }
320 ++leit;
321 }
322 return 0;
323 }
324
325 //! @brief this class only handle Ai+Bi=Ci
326 template <typename T>
327 class gauss_eliminator{
328 BOOST_MPL_ASSERT_MSG((boost::is_same<T,double>::value ) ||
329 (boost::is_same<T,float>::value ),
330 NON_FLOAT_TYPES_ARE_NOT_ALLOWED, (void));
331 public:
332 /**
333 * @brief construct gauss_eliminator class
334 *
335 * @param nodes input nodes
336 * @param node_flag input node_flag which will be tagged as true if the
337 * corresponding node is known
338 */
339 gauss_eliminator(zjucad::matrix::matrix<T> & nodes,
340 std::vector<bool> & node_flag)
341 :nodes_(nodes), node_flag_(node_flag){
342 idx2equation_.resize(nodes_.size());
343 }
344
345 /**
346 * @brief add equation to gauss_eliminator, every time an equation is added,
347 * eliminate function is called.
348 *
349 * @param input equation
350 * @return int
351 */
352 int add_equation(const equation<T> & e);
353
354 /**
355 * @brief This function will start to eliminate equations above all added equations
356 *
357 * @return int return 0 if works fine, or return non-zeros
358 */
359 int eliminate();
360
361
362 /**
363 * @brief update the equation, it will check all variant, if a variant is
364 * already known, update this equation
365 *
366 * @param eq input equation
367 * @return int return 0 if nothing changes, or retunr 1;
368 */
369 int update_equation(equation<T> & eq);
370
371 private:
372 zjucad::matrix::matrix<T> & nodes_;
373 std::vector<bool> & node_flag_;
374 std::list<equation<T> > es;
375
376 typedef typename std::list<equation<T> >::iterator equation_ptr;
377 std::vector<std::list<equation_ptr> > idx2equation_;
378
379 // this map store the smallest expression
380 typedef typename std::map<size_t, std::list<equation_ptr> >::iterator prime_eq_ptr;
381 std::map<size_t, std::list<equation_ptr> > prime_idx2equation_;
382 };
383
384 template <typename T>
385 int gauss_eliminator<T>::add_equation(const equation<T> & e){
386 es.push_back(e);
387 equation<T> & e_back = es.back();
388 for(typename equation<T>::eq_iterator eit = e_back.begin(); eit != e_back.end(); ){
389 if(node_flag_[eit->index]){
390 e_back.value() -= nodes_[eit->index] * eit->coefficient;
391 e_back.e_vec_.erase(eit++);
392 }else
393 ++eit;
394 }
395 if(e_back.state() == 0){// this equation is cleared
396 es.pop_back();
397 return 0;
398 }else if(e_back.state() == -1){
399 std::cerr << "# [error] strange conflict equation: " << std::endl;
400 std::cerr << e;
401 es.pop_back();
402 return __LINE__;
403 }
404 e_back.standardization();
405
406 for(typename equation<T>::eq_const_iterator it = e.begin();
407 it != e.end(); ++it){
408 equation_ptr end_ptr = es.end();
409 idx2equation_[it->index].push_back(--end_ptr);
410 }
411
412 equation_ptr end_ptr = es.end();
413 prime_idx2equation_[e_back.get_prim_idx()].push_back(--end_ptr);
414 eliminate();
415 return 0;
416 }
417
418 template <typename T>
419 int gauss_eliminator<T>::update_equation(equation<T> & eq)
420 {
421 for(typename std::list<expression<T> >::iterator it = eq.e_vec_.begin();
422 it != eq.e_vec_.end(); ){
423 const expression<T> & exp = *it;
424 if(node_flag_[exp.index]){
425 eq.value() -= exp.coefficient * nodes_[exp.index];
426 eq.e_vec_.erase(it++);
427 }else
428 ++it;
429 }
430 eq.standardization();
431 return 0;
432 }
433
434 template <typename T>
435 int gauss_eliminator<T>::eliminate()
436 {
437 std::cerr << std::endl;
438 while(1){
439 bool is_modified = false;
440
441 for(prime_eq_ptr ptr = prime_idx2equation_.begin();
442 ptr != prime_idx2equation_.end();)
443 {
444 std::list<equation_ptr> & dle = ptr->second;
445 if(dle.empty()) {
446 prime_idx2equation_.erase(ptr++);
447 is_modified = true;
448 continue;
449 }else if(dle.size() == 1){ // contain only one equation
450 const int state_ = dle.front()->state();
451 if(state_ == 0) { // cleared
452 es.erase(dle.front());
453 continue;
454 }else if(state_ == -1){ // conflict equation
455 std::cerr << "# [error] conflict equation " << std::endl;
456 return __LINE__;
457 }else if(state_ == 1){ // finial variant
458 const equation<T> & eq = *dle.front();
459 const T &value_ = eq.value();
460 // prime index's coefficient should be 1
461 assert(fabs(eq.e_vec_.front().coefficient - 1) < 1e-6);
462 const size_t index = eq.get_prim_idx();
463 if(node_flag_[index]){
464 if(fabs(nodes_[index] - value_) > 1e-6){
465 std::cerr << "# [error] conficts happen, node " << index
466 << " has different value " << nodes_[index] << ","
467 << value_ << std::endl;
468 return __LINE__;
469 }
470 }else{
471 // update corresponding equations with the new node value
472 nodes_[index] = value_;
473 node_flag_[index] = true;
474 std::list<equation_ptr> & node_linked_eq = idx2equation_[index];
475 for(typename std::list<equation_ptr>::iterator leqit =
476 node_linked_eq.begin(); leqit != node_linked_eq.end(); ++leqit){
477 equation<T> & eq = *(*leqit);
478 eq.update(index, nodes_[index]);
479 eq.standardization();
480 }
481 node_linked_eq.clear();
482 prime_idx2equation_.erase(ptr);
483 is_modified = true;
484 }
485 }
486 ++ptr;
487 }else{
488 assert(dle.size() > 1);
489 // this prime_index point to several equations,
490 // which sould be eliminated
491 typename std::list<equation_ptr>::iterator begin = dle.begin();
492 typename std::list<equation_ptr>::iterator first = begin++;
493 // to keep each prime index linked only one equation
494 for(typename std::list<equation_ptr>::iterator next = begin;
495 next != dle.end();){
496 // to eliminate the prime index, each equation minus the first one
497 *(*next) -= *(*first);
498 (*next)->standardization();
499 const size_t prim_index = (*next)->get_prim_idx();
500 assert(prim_index >= (*first)->get_prim_idx());
501 prime_idx2equation_[prim_index].push_back(*next);
502 dle.erase(next++);
503 }
504 ++ptr;
505 } // end else
506 }
507 if(!is_modified) break;
508 }
509 return 0;
510 }
511 }
512 }
513 #endif // GAUSS_ELIMINATION_H
下面有個簡單的測試例子：
//2x2+2x1+2x1=2
//x1=1
//5x1+x2=2
1 int main()
2 { typedef double val_type;
3 matrix<val_type> node = zeros<val_type>(4,1);
4 vector<bool> node_flag(4,false);
5
6 jtf::algorithm::gauss_eliminator<val_type> ge(node, node_flag);
7
8 {
9 equation<val_type> eq;
10 eq.add_expression(jtf::algorithm::make_expression(static_cast<size_t>(1),
11 static_cast<val_type>(2)));
12 eq.add_expression(jtf::algorithm::make_expression(static_cast<size_t>(0),
13 static_cast<val_type>(2)));
14 eq.add_expression(jtf::algorithm::make_expression(static_cast<size_t>(0),
15 static_cast<val_type>(2)));
16 eq.value() = 2;
17 ge.add_equation(eq);
18 }
19 {
20 equation<val_type> eq;
21 eq.add_expression(jtf::algorithm::make_expression(static_cast<size_t>(0),
22 static_cast<val_type>(1)));
23 eq.value() = 1;
24 ge.add_equation(eq);
25 }
26 {
27 equation<val_type> eq;
28 eq.add_expression(jtf::algorithm::make_expression(static_cast<size_t>(0),
29 static_cast<val_type>(5)));
30
31 eq.add_expression(jtf::algorithm::make_expression(static_cast<size_t>(1),
32 static_cast<val_type>(1)));
33
34 eq.value() = 2;
35 ge.add_equation(eq);
36 }
37 for(size_t t = 0; t < node_flag.size(); ++t){
38 if(node_flag[t] == true){
39 cerr << "# node " << t << " = " << node[t] << endl;
40 }else
41 cerr << "# node " << t << " unknown." << endl;
42 }
43 return 0;
44 }
這段code在g++4.6.3上編譯通過，如果有更加方便簡單的方法，請大家指教~
浙公網安備 33010602011771號