I'm in a situation where I have to be able to solve an unknown number of simultaneous equations, there could be from two variables to a near unlimited number (practical limits will probably be in the 10-15 area, but this program contains a lot of random events, so ther could potentially be much more).

The equations are all simple linear equations, and there will always be N equations, where N is the number of variables to solve. So the situation is not overly complex, however I need to do this fast (possibly millions of times in s single execution) so I want to try and find the fastest way possible of doing this, I've been pointed to a book called "Introduction to Algorithms", by T.H. cormen, C.E,Leisserson, R.L. Rivest.

And that's given by a psuedocode algorithm that is a factor of 3 faster than methods such as using a matrix inverse. The problem comes when I try to translate the code to C.

I realise this is a pretty complex topic, and most people probably won;t even want to bother wading through the code below, but I'm sorta out of options :(. If anyone's got some suggestions on other approaches or good sources of relevant information, I'd REALLY love to hear about it.

Anyway, this is the psuedo-code I have:

if anyone actually has the book, it's a fair bit easier to understand there, with the characters and formatting I just can't get in a forum post.Code:`LUP-solve (L,U,Pi,b)`

n = rows[L]

for i = 1 to n

do y[i] = b[Pi[i]] - sum(j=1 to i-1) of L[i][j]y[i]

for i = n downto 1

do x[i] = (y[i]-sum(j=i+1 to n) of u[i][j]x[j])/u[i][i]

return x

LUP-Decomposition(A)

n = rows[A]

for i = 1 to n

do Pi[i] = i

for k = 1 to n-1

do p = 0

for i = k to n

do if |A[i][k]| > p

then p=|A[i][k]|

k' = i

if p=0

then error "singular matrix"

exchange Pi[k] <> Pi[k']

for i=1 to n

do exchange A[k][i] <> A[k'][i]

for i=k+1 to n

do A[i][k] = A[i][k]/A[k][k]

for j =k+1 to n

do A[i][j] = A[i][j] - A[i][k]A[k][j]

My attempt at this code is:

I've compressed the two functions into one here, though that's not a problem, in the initial version, the decomposition runs, then solve, I've done the same thing, just in one function call.Code:

int LUPsolve(int rows,float A[rows][rows],float b[rows],float X[rows])

{

int P[rows];

int n,m,p,q,i,j;

int k=0,kd=0;

int T=0;

float t=0;

float Y[rows];

//decompose

for(i=1;i<rows;i++)

{

P[i]=i;

}

for(k=1;k<rows-1;k++)

{

p=0;

for(i=k;i<rows;i++)

{

t=A[i][k];

if(A[i][k]<0)t=-1*t;

if(t>p)

{

p=t;

kd=i;

}

}

if(p==0)

{

printf("singular Matrix, no possible solutions\n");

return 1;

}

T=P[kd];

P[kd]=P[k];

P[k]=T;

for(i=1;i<rows;i++)

{

t=A[kd][i];

A[kd][i]=A[k][i];

A[k][i]=t;

}

for(i=k+1;i<rows;i++)

{

A[i][k]=A[i][k]/A[k][k];

for(j=k+1;j<rows;j++)

{

A[i][j]=A[i][j]-A[i][k]*A[k][j];

}

}

}

//now solve

for(n=1;n<rows;n++)

{

t=0;

for(m=1;m<n-1;m++)

{

t+=A[n][m]*Y[m];

}

Y[n]=b[P[n]]-t;

}

for(n=rows-1;n>=1;n--)

{

t=0;

for(m=n+1;m<rows;m++)

{

t+=A[n][m]*X[m];

}

X[n]=(Y[n]-t)/A[n][n];

}

//X now contains the solution;

return 0;

}

Also, though the psuedo-code for LUP-solve doesn't show it (decom does) the matrices L and U are both combined into A, so that is not the problem with the code. The indexing from 1 rather than 0 may be a problem, but when I go through and edit it to index from 0, it doesn;t work either, so I just keep it like this to maintain continuity with the book.

the input I'm tesing it with is:

a={{0,0,0},{0,4,2},{0,3,-2}}

b={0,6,8}

rows=3

done on paper, the output should be 2 and -11, but I actually get 2.64 & -2.28

anyone got an idea what's wrong/any pointers to more resources/pre-coded functions/anything remotely useful. thisis something I really want to be able to do, but my own brain's pretty much worn out on this :(

big thanks to anyone who's still reading by now :)

-mark