X-Git-Url: http://git.tdb.fi/?a=blobdiff_plain;f=source%2Flinal%2Fmatrixops.h;h=0988accd6d44798ab906851fc180acf4455cba3e;hb=ff8619c1c440fd1a7c8deaca763a33b8e1d53b5e;hp=790dfbfaf694724f6040555077b56e920ff0c4f3;hpb=5797fe0a5296952cb8f8643fdc6cabddee19a554;p=libs%2Fmath.git diff --git a/source/linal/matrixops.h b/source/linal/matrixops.h index 790dfbf..0988acc 100644 --- a/source/linal/matrixops.h +++ b/source/linal/matrixops.h @@ -3,6 +3,7 @@ #include #include +#include namespace Msp { namespace LinAl { @@ -16,44 +17,44 @@ public: template -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; iabs(m.element(pivot, i))) + for(unsigned j=i+1; jabs(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; j0; ) + 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