1 #ifndef MSP_LINAL_MATRIXOPS_H_
2 #define MSP_LINAL_MATRIXOPS_H_
10 class not_invertible: public std::domain_error
13 not_invertible(): domain_error(std::string()) { }
14 virtual ~not_invertible() throw() { }
19 inline T &invert_matrix(T &m, T &r)
21 typedef typename T::ElementType V;
24 for(unsigned i=0; i<m.rows(); ++i)
27 for(unsigned j=i+1; j<m.rows(); ++j)
28 if(abs(m.element(j, i))>abs(m.element(pivot, i)))
31 if(m.element(pivot, i)==V(0))
32 throw not_invertible();
36 m.exchange_rows(i, pivot);
37 r.exchange_rows(i, pivot);
40 for(unsigned j=i+1; j<m.rows(); ++j)
42 V a = -m.element(j, i)/m.element(i, i);
47 V a = V(1)/m.element(i, i);
52 for(unsigned i=m.rows(); i-->0; )
53 for(unsigned j=i; j-->0; )
54 r.add_row(i, j, -m.element(j, i));