faif
GaussEliminator.h
1 #ifndef FAIF_GAUSS_ELIMINATOR
2 #define FAIF_GAUSS_ELIMINATOR
3 
4 // plik zawiera deklaracje przeksztalcenia - prosty eliminator Gaussowski (niestety O(N^3) )
5 
6 #if defined(_MSC_VER) && (_MSC_VER >= 1400)
7 //msvc8.0 generuje warning dla boost::numeric::matrix
8 #pragma warning(disable:4996)
9 #endif
10 
11 #include <cassert>
12 #include <boost/numeric/ublas/vector.hpp>
13 #include <boost/numeric/ublas/matrix.hpp>
14 
15 namespace faif {
16 
17  template<typename T>
18  //parametry przekazywane przez referencje, sa niszczone podczas obliczen
19  boost::numeric::ublas::vector<T> GaussEliminatorRef(boost::numeric::ublas::matrix<T>& m,
20  boost::numeric::ublas::vector<T>& y) {
21 
22  // std::assert( y.size() == m.size1() && y.size() == m.size2() );
23 
24  int n = static_cast<int>(y.size() ); //wielkosc
25 
26  //tworzy macierz trojkatna
27  for (int i = 0; i < n; ++i) {
28  int max = i; //maksymalny element
29  for (int j = i + 1; j < n; ++j)
30  if( m(j,i) > m(max,i) )
31  max = j;
32 
33  for (int j = 0; j < n; ++j)
34  std::swap( m(max,j), m(i,j) );
35  std::swap( y(max), y(i) );
36 
37  for(int k = i + 1; k < n; ++k)
38  y(k) -= m(k,i)/m(i,i) * y(i);
39  for (int j = n-1; j >= i; --j)
40  for (int k = i + 1; k < n; ++k)
41  m(k,j) -= m(k,i)/m(i,i) * m(i,j);
42  }
43 
44  //eliminacja - oblicza wyniki
45  for (int i = n - 1; i >= 0; --i) {
46  y(i) = y(i)/m(i,i);
47  m(i,i) = 1;
48  for (int j = i - 1; j >= 0; --j) {
49  y(j) -= m(j,i)*y(i);
50  m(j,i) = 0;
51  }
52  }
53  return y;
54  }
55 
56  template<typename T>
57  //kopiuje parametry, wywoluje GaussEliminatorRef, parametry nie sa niszczone
58  boost::numeric::ublas::vector<T> GaussEliminator(const boost::numeric::ublas::matrix<T>& m,
59  const boost::numeric::ublas::vector<T>& y) {
60  typename boost::numeric::ublas::matrix<T> mm(m);
61  typename boost::numeric::ublas::vector<T> yy(y);
62  return GaussEliminatorRef(mm,yy);
63  }
64 
65 
66 } //namespace faif
67 
68 #endif //FAIF_GAUSS_ELIMINATOR
Definition: Chain.h:17