## Friday, December 11, 2009

### Gauss–Jordan Elimination in C++

This is my first attempt to solve linear systems with unique solutions in O(N3).

`Algorithm:1. Read augmented matrix for the system.2. Perform Gaussian Elimination to get the matrix in row echelon form.3. Transform the matrix to reduced row echelon form.4. Write results.`

Preview:

[click for a clear view]

My C++ solution (change MAX to desired size, it is the maximum limit of how many equations it can handle, change the macro DEBUG to "if(1)" if you want to see the intermediate steps):
`#include <iostream>#include <cstdlib>#include <cassert>#include <algorithm>using namespace std;#define DEBUG if(0)#define MAX 5struct MAT {    int R[MAX+1];} M[MAX];int N;int gcd(int a, int b) {    while(b) b ^= a ^= b ^= a %= b;    return a;}bool comp(MAT A, MAT B) {    for(int k=0; k<N; k++) {        if(A.R[k] > B.R[k]) return true;        if(A.R[k] < B.R[k]) return false;    }    return true;}void print(void) {    int i, j;    cout << "..." << endl;    for(i=0; i<N; i++) {        for(j=0; j<N; j++) cout << M[i].R[j] << '\t';        cout << M[i].R[j] << endl;    }    cout << "..." << endl;}void modify(MAT &A, MAT B, int a, int b) {    for(int r=0; r<=N; r++) A.R[r] = a*B.R[r] - b*A.R[r];}void eliminate(void) {    int i, g, k;    for(i=0; i<N-1; i++) {        if(!M[i].R[i]) {            sort(&M[i], M+N, comp);            assert(M[i].R[i]);        }        for(k=i+1; k<N; k++) {            g = gcd(abs(M[i].R[i]), abs(M[k].R[i]));            modify(M[k], M[i], M[k].R[i]/g, M[i].R[i]/g);        }        DEBUG print();    }}void reduce(void) {    int i, g, k;    for(i=N-1; i; i--) {        for(k=i-1; k>=0; k--) {            g = gcd(abs(M[i].R[i]), abs(M[k].R[i]));            modify(M[k], M[i], M[k].R[i]/g, M[i].R[i]/g);        }        DEBUG print();    }}void printsol(void) {    double h, l;    cout << "Solve for " << N << " variables:" << endl;    for(int i=0; i<N; i++) {        assert(M[i].R[i]);        h = M[i].R[i];        l = M[i].R[N];        cout << 'x' << i+1 << " = " << l/h << endl;    }    cout << "..." << endl;}int main(void) {    int i, j, g;    while(cin >> N) {        if(!N) break;        memset(M, 0, sizeof(M));        for(i=0; i<N; i++) {            for(j=0; j<N; j++) cin >> M[i].R[j];            cin >> M[i].R[j];            for(j=g=0; j<=N; j++) g = gcd(abs(M[i].R[j]), g);            for(j=0; j<=N; j++) M[i].R[j] /= g;        }        eliminate();        reduce();        printsol();    }    return 0;}`

However, I'm trying to do it in a more efficient way. I'll post that when I'm done.

1. WOWOWOWOWOWOWOW!

2. Vaiya Great...

3. nah, I failed to find any better way other than some nasty coding optimisation...
any one knows a better way?

4. Yeah, nor do I.

5. Plz provide some links of problems related to this topic..

1. Will try to when I have some time :)

6. How is this O(N^2) ? Isn't this O(N^3) ?

1. Of course its N^3, must be a typo. Ill fix that.

2. Fixed. Thanks :)