- for(unsigned i=0; i<S; ++i)
- {
- unsigned pivot = i;
- for(unsigned j=i+1; j<S; ++j)
- if(abs(this->element(j, i))>abs(this->element(pivot, i)))
- pivot = j;
-
- if(this->element(pivot, i)==T(0))
- throw not_invertible();
-
- if(pivot!=i)
- {
- this->exchange_rows(i, pivot);
- r.exchange_rows(i, pivot);
- }
-
- for(unsigned j=i+1; j<S; ++j)
- {
- T a = -this->element(j, i)/this->element(i, i);
- this->add_row(i, j, a);
- r.add_row(i, j, a);
- }
-
- T a = T(1)/this->element(i, i);
- this->multiply_row(i, a);
- r.multiply_row(i, a);
- }
-
- for(unsigned i=S; i-->0; )
- for(unsigned j=i; j-->0; )
- r.add_row(i, j, -this->element(j, i));
-