#include <cmath>
#include <stdexcept>
+#include <string>
namespace Msp {
namespace LinAl {
template<typename T>
-inline T &invert_matrix(T &m, T &r)
+inline T &gauss_jordan(T &m, T &r)
{
typedef typename T::ElementType V;
using std::abs;
- for(unsigned i=0; i<m.rows(); ++i)
+ for(unsigned i=0; i<m.columns(); ++i)
{
unsigned pivot = i;
- for(unsigned j=i+1; j<m.rows(); ++j)
- if(abs(m.element(j, i))>abs(m.element(pivot, i)))
+ for(unsigned j=i+1; j<m.columns(); ++j)
+ if(abs(m.element(i, j))>abs(m.element(i, pivot)))
pivot = j;
- if(m.element(pivot, i)==V(0))
+ if(m.element(i, pivot)==V(0))
throw not_invertible();
if(pivot!=i)
{
- m.exchange_rows(i, pivot);
- r.exchange_rows(i, pivot);
+ m.exchange_columns(i, pivot);
+ r.exchange_columns(i, pivot);
}
- for(unsigned j=i+1; j<m.rows(); ++j)
+ for(unsigned j=i+1; j<m.columns(); ++j)
{
- V a = -m.element(j, i)/m.element(i, i);
- m.add_row(i, j, a);
- r.add_row(i, j, a);
+ V a = -m.element(i, j)/m.element(i, i);
+ m.add_column(i, j, a);
+ r.add_column(i, j, a);
}
V a = V(1)/m.element(i, i);
- m.multiply_row(i, a);
- r.multiply_row(i, a);
+ m.multiply_column(i, a);
+ r.multiply_column(i, a);
}
- for(unsigned i=m.rows(); i-->0; )
+ for(unsigned i=m.columns(); i-->0; )
for(unsigned j=i; j-->0; )
- r.add_row(i, j, -m.element(j, i));
+ r.add_column(i, j, -m.element(i, j));
- return m = r;
+ return r;
}
} // namespace LinAl