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.rows(); ++i)
28 for(unsigned j=i+1; j<m.rows(); ++j)
29 if(abs(m.element(j, i))>abs(m.element(pivot, i)))
32 if(m.element(pivot, i)==V(0))
33 throw not_invertible();
37 m.exchange_rows(i, pivot);
38 r.exchange_rows(i, pivot);
41 for(unsigned j=i+1; j<m.rows(); ++j)
43 V a = -m.element(j, i)/m.element(i, i);
48 V a = V(1)/m.element(i, i);
53 for(unsigned i=m.rows(); i-->0; )
54 for(unsigned j=i; j-->0; )
55 r.add_row(i, j, -m.element(j, i));