1 #ifndef MSP_LINAL_MATRIXOPS_H_
2 #define MSP_LINAL_MATRIXOPS_H_
11 class not_invertible: public std::domain_error
14 not_invertible(): domain_error(std::string()) { }
15 virtual ~not_invertible() throw() { }
20 inline T &gauss_jordan(T &m, T &r)
22 typedef typename T::ElementType V;
25 for(unsigned i=0; i<m.columns(); ++i)
28 for(unsigned j=i+1; j<m.columns(); ++j)
29 if(abs(m.element(i, j))>abs(m.element(i, pivot)))
32 if(m.element(i, pivot)==V(0))
33 throw not_invertible();
37 m.exchange_columns(i, pivot);
38 r.exchange_columns(i, pivot);
41 for(unsigned j=i+1; j<m.columns(); ++j)
43 V a = -m.element(i, j)/m.element(i, i);
44 m.add_column(i, j, a);
45 r.add_column(i, j, a);
48 V a = V(1)/m.element(i, i);
49 m.multiply_column(i, a);
50 r.multiply_column(i, a);
53 for(unsigned i=m.columns(); i-->0; )
54 for(unsigned j=i; j-->0; )
55 r.add_column(i, j, -m.element(i, j));