LUP Decomposition (simultaneous equations)

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:

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]

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.

My attempt at this code is:

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;
}

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.

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

I don't know if this will be of any use, but here is some C++ code that might be helpful to look at.

Code:

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

This compiles? Is this C++? C99?

it compiles and runs in dev_c++ in C mode.

I know using "rows" all over th plpace there is not really a great way of doing things, but if I leave them out I get compile errors, something to do with variable arrays.

Big thanks for that link, lines 85-128 seem exactly what I'm looking for, shame I don;tunderstand all the differences between C/C++ but I should be able to make out enough of it to dig around for the problems. I've bookmarked that and I'll deal with it in the morenig when I'm a little more awake (11:20 Pm is not the best time for hits kinda stuff...).

thanks again

-mark

Mark, can you explain what U is. In the pseudo-code:
LUP-solve (L,U,Pi,b)

Now it appears that LUP-decomposition is done first, then LUP-solve. But in the pseudocode, U is never calculated. I'm assuming it's not the b array
b=[0,6,8]
as this seems to be passed last.

I know you said L and U are combined into A, but you're using A as if it's L in the pseudocode.

the pseudocode comes with 11 pages of further information detailing what's actually going on, a bit too much to repeat here, but I'll try and summarize what I can, though I should warn that I'm terrible at explanations:

L and U are triangular matrices, (sort of only half a matrix) for example:
L=
1 0 0
0.6 1 0
0.2 0.5 1

U=
5 6 3
0 1.4 2.2
0 0 -1.8

notice that they both have lots of zeros in the top right/bottom left corner.

the LUP decomposition returns the matrix A, which contains both of these matrices:

5 6 3
0.6 1.4 2.2
0.2 0.5 -1.8

both L and U can be reconstructed from A, as the only cells that are not passed directly are ones.

if you look at the sum statements in the LUP-Solve function, each sum is only summing one part of A.

for i = 1 to n
sum(j=1 to i-1)

implies that (for n=3) we sum cells

i=1: no cells
i=2: j=1
i=3: j=1, j=2

and for the second:

i=3:no cells
i=2:j=3
i=1:j=2,j=3

so if we write the matri, showing which sum read each cell:

_22
1_2
11_

they're both reading triangles from the same matrix and I think I just saw a potential problem...the diagonal in the middle isna't being read at all, even though it's got values in. I think I'll take a closer look at that, thanks for making me step through it in detail again

-mark

I'm really getting sick of silly errors, managed to have one in here that's been there for weeks.... at least I've got it now, the correct code is:

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];
  //initialise X,Y
  for (n=1;n<rows;n++)
    {
    X[n]=0;
    Y[n]=0;
    }
  //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(t<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;
}

the error was in the third last for statement, changing

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

to

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

now works fine, dunno how I missed that in the who-knows-how many times I worked through the code & psuedocode... silly errors seems to be a really big thing for me

-mark

Cool program. I ran your code with the values above and got 2 and -1. That was a tough bug to spot.