From: Mikko Rasa Date: Sun, 13 Mar 2022 19:32:47 +0000 (+0200) Subject: Make gauss_jordan operate on columns instead of rows X-Git-Url: http://git.tdb.fi/?a=commitdiff_plain;h=ff8619c1c440fd1a7c8deaca763a33b8e1d53b5e;p=libs%2Fmath.git Make gauss_jordan operate on columns instead of rows Since matrices are stored in column-major order, this makes the algorithm more amenable to possible simdification. --- diff --git a/source/linal/dynamicmatrix.h b/source/linal/dynamicmatrix.h index 1b191a1..de4511c 100644 --- a/source/linal/dynamicmatrix.h +++ b/source/linal/dynamicmatrix.h @@ -45,9 +45,9 @@ public: DynamicMatrix &operator+=(const DynamicMatrix &); DynamicMatrix &operator-=(const DynamicMatrix &); - DynamicMatrix &exchange_rows(unsigned, unsigned); - DynamicMatrix &multiply_row(unsigned, T); - DynamicMatrix &add_row(unsigned, unsigned, T); + DynamicMatrix &exchange_columns(unsigned, unsigned); + DynamicMatrix &multiply_column(unsigned, T); + DynamicMatrix &add_column(unsigned, unsigned, T); DynamicMatrix &invert(); }; @@ -252,38 +252,38 @@ inline bool operator==(const DynamicMatrix &m1, const DynamicMatrix &m2) } template -inline DynamicMatrix &DynamicMatrix::exchange_rows(unsigned i, unsigned j) +inline DynamicMatrix &DynamicMatrix::exchange_columns(unsigned i, unsigned j) { - if(i>=rows_ || j>=rows_) - throw std::out_of_range("DynamicMatrix::exchange_rows"); + if(i>=columns_ || j>=columns_) + throw std::out_of_range("DynamicMatrix::exchange_columns"); using std::swap; for(unsigned k=0; k -inline DynamicMatrix &DynamicMatrix::multiply_row(unsigned i, T s) +inline DynamicMatrix &DynamicMatrix::multiply_column(unsigned i, T s) { - if(i>=rows_) - throw std::out_of_range("DynamicMatrix::multiply_row"); + if(i>=columns_) + throw std::out_of_range("DynamicMatrix::multiply_column"); for(unsigned k=0; k -inline DynamicMatrix &DynamicMatrix::add_row(unsigned i, unsigned j, T s) +inline DynamicMatrix &DynamicMatrix::add_column(unsigned i, unsigned j, T s) { - if(i>=rows_ || j>=rows_) - throw std::out_of_range("DynamicMatrix::exchange_rows"); + if(i>=columns_ || j>=columns_) + throw std::out_of_range("DynamicMatrix::exchange_columns"); for(unsigned k=0; k @@ -227,27 +227,27 @@ inline bool operator==(const Matrix &a, const Matrix &b) } template -inline Matrix &Matrix::exchange_rows(unsigned i, unsigned j) +inline Matrix &Matrix::exchange_columns(unsigned i, unsigned j) { using std::swap; - for(unsigned k=0; k -inline Matrix &Matrix::multiply_row(unsigned i, T s) +inline Matrix &Matrix::multiply_column(unsigned i, T s) { - for(unsigned k=0; k -inline Matrix &Matrix::add_row(unsigned i, unsigned j, T s) +inline Matrix &Matrix::add_column(unsigned i, unsigned j, T s) { - for(unsigned k=0; kabs(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 r; }